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Chapter 1 
Introduction 



Quantum many-body problems in condensed matter physics are a context 
of everlasting interest and relentless investigation in physical research. The 
macroscopic amount of interacting degrees of freedom is such that even the 
simplest models become extremely hard problems, in exact description as 
well as in pcrturbative regime. The analytical and computational complex- 
ity of many-body physics is deeply rooted in the foundations of quantum 
mechanics themselves: the hardness of these problems undergoes a scaling- 
law with the number of elementary constituents (size) of the system L, which 
is typically much more abrupt than extensive behavior; it usually grows expo- 
nentially with L. In physical literature a variety of models which have proven 
to be particularly suitable for analytical study was developed, e.g. due to 
some peculiar local structure or to some wide symmetry-group invariance; 
however, the large majority of known non-perturbative Hamiltonians mani- 
fest no attitude towards analytical simplification and must be faced head-on 
with numerical techniques. Efficient simulation methods for condensed mat- 
ter settings are many, often capable of integrating any bit of theoretical 
knowledge, then crossing the remaining gap with computation. At the same 
time, when no previous hint from theory is available, the exact problem turns 
hard again, and computational costs scale fast with system specifics, so that 
only very small sizes L are manageable for practical purposes. 

The family of variational algorithms has always been regarded as one 
of the most natural and promising paths in order to address many-body 
problems at zero temperature: as ground states of Hamiltonians are minima 
of the spectrum, searching them through variational procedures seems most 
appealing. Yet the crucial point of any variational paradigm is the capability 
of reducing the whole amount of degrees of freedom into a small number 
of effective, important ones: these must embed all the relevant physics of 
the target state, while requiring a limited number of numerical resources. 
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Such primary descriptors, or variational parameters (also known as reaction 
coordinates in some contexts) need to be identified and discriminated from 
the non-influential ones. This is clearly a delicate issue, especially if no 
knowledge on the model is available. In other words, the selection a priori of 
some appropriate basis of variational wavefunctions is the fundamental step 
to undertake, determining the faithfulness and efficiency (and thus overall 
success) of any variational algorithm we might want to develop. 

In this thesis, we will focus on a very general family of variational wave- 
functions, whose main peculiarity is that their descriptors/parameters are 
tailored according to simple linear algebraic relations. The computational 
power and success of these tools descends from arguments that were born 
within quantum information framework: entanglement [1]. Quantum entan- 
glement is indeed a resource, but it is also a measure of internal correlations 
in multipartite systems. Once we characterized general entanglement prop- 
erties of many-body ground states, then by controlling entanglement of a 
variational trial wavefunction we can exclusively address physical states, and 
disregard non-physical states, even before the simulation takes place. This 
is the central concept which Tensor Network architectures are based upon. 

Historically, the realization and profound understanding of this class of 
states, was possible only after the introduction of Density Matrix Renormal- 
ization Group (DMRG), curiously, an algorithm which is not formulated in 
variational terms at all. 

1.1 The age of Density Matrix 
Renormalization Group 

The idea of adapting a Renormalization Group algorithm to a lattice Den- 
sity Matrix was proposed by Steven R. White, considered undoubtedly one 
of the founders of the DMRG methods. In his first approaches [2] to the 
technique, he was inspired by a paper of K.W. Wilson |3] where a numerical 
renormalization group paradigm is applied to the Hamiltonian of a Kondo 
problem. 

The simple, yet brilliant, idea behind White's formulation of DMRG was 
to replace the traditional procedure of renormalization group, which acted 
extensively in the real space, and thus actually performing a coarse-graining 
transformation upon the system, with a scheme that applied extensively in 
the Hilbert space dimension itself, leading to a site-by-site renormalization 
scheme. In practice, assume that we are describing the state of a given 
portion of the system (in terms of a density matrix) with a fixed amount 
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of computational resources D. Now when another single constituent (a site) 
is added to the subsystem, the resulting dimension grows linearly with the 
local dimension d of the new component, ~ dD. Then renormalization is 
performed, allowing us to represent the new subsystem with the same initial 
amount of resources D, cutting the least relevant density matrix-eigenstates 
out of the description. Obviously, such operations still manifest a group 
structure, and they are summoned every time the density matrix dimension 
(and not the real-space size) increases of a given factor d: thus DMRG. 

The great amount of credit and interest gathered by DMRG is surely 
due to its outstanding successfulness for low-dimensionality quantum sys- 
tems. In particular, for one- dimensional (open-boundary) systems, DMRG 
achieved variational precisions (compared to experiment, and theory when- 
ever possible) that challenged other simulation approaches. At the end of 
the '90s, it was considered probably the most powerful numerical method to 
address ID problems, with practically no knowledge on the model required a 
priori. In literature, DMRG picture has been exploited in several settings of 
both physics and quantum chemistry, and numerous variants of its original 
formulation were proposed [H |5]; in the end the basic idea was proven to be 
winning, even though within its dimensionality limits. 

1.2 The advent of Matrix Product States 

The DMRG concept was quite renown, but it was in 2004 that the in-depth 
reason of its success was fully understood: when F. Verstraete, D. Porras and 
J. I. Cirac started to investigate quantum states built via a DMRG algorithm 
under a quantum information perspective [6J. In fact, they realized that 
DMRG states had a strict equivalence relation with finitely correlated states, 
i.e. lattice states whose entanglement is upper-bounded by an arbitrary finite 
value, which does not scale with the system size. 

Moreover, ground states of short-ranged (non-critical) Hamiltonians, have 
been known for quite some time to satisfy the so-called area-law of entan- 
glement [3 [HI [9]. This general rule was developed in quantum information 
contexts, but carries important physical prescriptions. It claims that the 
partition entanglement of a non-critical ground state does not scale as the 
volume of the parted spatial region ~ L*^ (which is the typical entanglement 
scaling for random states), but rather with the parting surface ~ 
where #iV is the number of spatial dimensions. 

It is clear that for ID systems, the correct area law is given by a non- 
scaling constant ~ L^: indeed entanglement of ID non-critical ground states 
typically saturates to a finite bound. This means, in turn, that ground states 
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of ID non-critical Hamiltonians are finitely-correlated states, and therefore, 
that DMRG procedure can approximate them with arbitrary precision, and 
their entanglement as well. 

Moreover, ref. ^ has shown that finitely-correlated states allow a direct, 
simple, immediate algebraic representation in terms of a product of matrices, 
each of these matrices storing all the information related to a single renor- 
malization process. This provides a one-to-one local correspondence between 
DMRG, and these Matrix Product State (MPS) [lOl [IH US US] which, as 
they allow an explicit analytic expression, actually form a class of tailored 
variational wavef unctions. 

Although not-homogeneously perceived by the condensed-matter physics 
community, this discovery was definitely a breakthrough, for several rea- 
sons. First, a variational formulation of DMRG opened the possibility for 
new numerical strategies based on finitely-correlated states, so that several 
minimum- search algorithms could be applied, but still using just the right 
amount of necessary computational resources. Secondly, the algebraic MPS 
expression provided faster ways to access physical information, and, at the 
same time, it allowed innovative problem-solving options even from the an- 
alytical point of view jT3]. Finally, it is easy to generalize the MPS concept 
to suit physical settings other than ID, and still taking care of appropriate 
area-laws. This argument lead, for instance, to the design of Product of 
Entangled Pairs States (PEPS) [161 [15]. 

1.3 Entanglement Renormalization 

The understanding of the relationship between DMRG, MPS, and finitely- 
correlated states provided an unquestionable paradigm for dealing numeri- 
cally with non-critical ID system in a fully-contextualized theoretical frame- 
work. Despite the fair success of adapting these algorithms to critical prob- 
lems as well (although in a hand- waving and unnatural way), it was guessed 
that due to the presence of scale-invariance symmetry, an old-fashioned real- 
space renormalization group would be more appropriate to simulate strongly 
correlated systems. Indeed, a state which is locally stationary under the 
action of a coarse-graining transformation would definitely be scale invari- 
ant. At the same time Wilson's numerical RG had the unpractical feature 
of suffering loss of short-range detailed structure, identified by translational 
instability of entanglement. 

An intriguing proposal to work around this trouble was introduced by 
G. Vidal in 2007, who fist spoke about Entanglement Renormalization [T7] . 
The idea is intuitive, yet very effective: we are still performing a real-space 
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renormalization group, but prior to the renormalization process itself, we 
apply a quasi-local unitary transformation, whose purpose is to decrease the 
correlation among regions who are to be renormalized separately. Since the 
goal of these unitary gates is to absorb, and thus store, entanglement out of 
the pre-RG state, they are commonly called disentanglers. The disentangling 
operation is then scheduled before every real-space RG operation takes place; 
in the end, such entanglement renormalization acts at every lenghtscale, since 
every RG step performs actually a scale transformation. 

Similarly to DMRG and MPS, these entanglement-RG states have a vari- 
ational counterpart as well. Precisely, it is possible to define a class of tai- 
lored variational wavefunctions, whose descriptors are tied together by linear 
relations, reproducing exactly the entire set of those states. Such states 
are thus called Multiscale Entanglement Renormalization Ansatz (MERA) 
[T8l flQl l20l [21] , and manifest a natural attitude towards describing strong 
correlation and criticality [22j. 

MPS, MERA, PEPS, are different classes of variational states sharing 
some important attributes: they are able to capture interesting physics, 
yet they require a small, manageable number of parameters, with simple 
algebraic rules and direct access to relevant physical information. Physi- 
cists started to regard them as belonging together to a larger, comprehensive 
family of tailored variational states, whose entanglement can be directly con- 
trolled through the selection of a related graph geometry. This is the concept 
of Tensor Network states [2S1 123 ESI I2S] ■ 

1.4 Outline 

The thesis is organized as follows: 

• In chapter [2] we will present an in-depth review on Matrix Product 
States, in ID open-boundary conditions settings. We will show how 
the MPS analytical expression is derived by the DMRG algorithm, 
show its entanglement bounds and sketch their algebraic manipulation 
features. We will explain how to achieve physical information on these 
states in a computationally-fast scheme, and present some protocols to 
simulate ground states. 

• In chapter [3] we will generalize the concept of MPS to periodic bound- 
ary conditions systems, and discuss how translational homogeneity of 
the representation allows us to well-define the thermodynamical limit 
for MPS. This will be the proper setting to show that MPS manifest 
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natural non-criticality, whose signature is an exponential decay of two- 
point correlations. We will then present possible generalizations to the 
PBC case of MPS minimization algorithms, focusing on some tricks of 
the trade useful to speed-up and stabilize the procedure. 

• In chapter |4] we will explain how the MPS architecture can be gener- 
alized to more complex Tensor Network geometries. We will investi- 
gate the entanglement properties of Tensor Network states, and some 
of their common algebraic features, like efficient contraction schemes 
or the adaptability to fermionic contexts. We will also present some 
remarkable subclasses of Tensor Networks, like PEPS and CPS, and 
discuss on how they relate to each other. 

• In chapter [5] we will put our attention on Trees and MERA, two classes 
of Tensor Networks that share most of their main features. We will 
show their natural attitude to describe critical states in ID, identified 
both by a logarithmic violation of the area law, and more importantly 
by manifesting power-law decaying correlations. Critical exponents, as 
well as the TTN/MERA state properties in the thermodynamical limit, 
are completely characterizable by adopting a completely positive trace- 
preserving map formalism. We will investigate other general properties 
of such TN-architectures, e.g. the possibility to construct a parent 
Hamiltonian. 



1.4.1 Original content 

Here I will list the original contribution I developed personally, either as 
brand-new material or as a reinterpretation of previous knowledge, during 
my Philosophiae Doctorateship. 



Section 2.11 analytical MPS representation of Slater Determinants, 



many-body basis change, and configuration interaction states. 



Section 4.8.1 matrix product representation for correlator product 



states (Jastrow factors). 

Most analytical results of chapter [5] scaling properties for TTN, parent 
Hamiltonians, fiuctuations, boundaries, hybrid geometries; and respec- 
tive generalizations to MERA. 



Part of this research was published in [271 [281 129] . 



Chapter 2 

Matrix Product States 



It was in 2004 that computational physicists started to consider Density Ma- 
trix Renormahzation Group according to a Quantum Information perspec- 
tive [0] ; they realized that it is possible to understand DMRG in a variational 
sense, in which the role played by entanglement and quantum correlation is 
clear. Indeed, as a quantum many-body state achieved by DMRG proce- 
dure is uniquely defined by the renormahzation transformations (intended as 
endomorphisms upon the density matrices space) one can regard such trans- 
formations as variational elements, and every single choice of those elements 
defines a state within a set of tailored variational wavef unctions. Nicely 
enough, it was discovered that such states allow a simple and immediate 
analytical description, where their many-body wavefunction, wrote upon a 
product basis of one-body levels, appears just as a product of variational 
matrices, thus leading to the name of Matrix Product States (MPS). As a 
matter of fact, such transparent description allowed research to further in- 
vestigate the properties of these states, leading to a deeper understanding of 
DMRG as well, and in the end the new knowledge served well the purpose of 
gaining more computational power in simulations, through a wider range of 
algebraic manipulations and the adaptability of variational-based algorithms. 

In the end, it was the Matrix Product State picture that helped to un- 
derstand the deep physical reason of DMRG successfulness in ID. Indeed, 
MPS have proven to be in tight relation with ID finitely correlated states 
[301 EI], and in turn this set is known to include ground states of short- 
range interacting non-critical Hamiltonians. Such argument holds not only 
for finite, isolated systems, but extends naturally to open and/or thermody- 
namical limit systems (as long as a zero temperature can be defined), allowing 
DMRG/MPS to succeed even in these cases. 

It is important to point out that Matrix Product State methods can be 
successfully adopted for dealing with fermionic systems, and they naturally 
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avoid the sign problem which, in stead, is a major issue on other variational 
fermionic algorithms like Montecarlo. Finally, MPS are not merely a numer- 
ical tool, they have proven in various contexts to be fundamental to address 
analytically several condensed matter models [321 [T3] . 

2.1 Matrix Product State construction 
from DMRG 

Following the formalism of |2], [10], we start with a one-dimensional lattice, 
L (length) being the total number of sites, d (local Hilbert dimension) the 
number of levels per site, and where open boundary conditions (OBC) are 
chosen for simplicity. Let us assume that we are describing a given quantum 
state of this system obtained via a DMRG algorithm: D <^ is the maximal 
number of states allowed for the each renormalization. Now let i < L be 
the site where the last density matrix renormalization was applied in the 
algorithm (while moving, say, to the right), this means that we know the 
reduced density matrix of the state p^, involving all the sites from i to the 
leftmost, but we only have access to its renormalized form: namely in stead of 
keeping in memory all its d^ eigenvectors \Lj)f and their relative probability 

Pj Pj = 1) in decreasing order pj > Pj+i), only the D < d^ oi such 
vectors are kept, of course those with highest probabihty: 

D 

Pi 



-J2p^\^^^^^^^\^ ^^^"^^ = (2-1) 



k=lPk 

ensures that the new statistic pj is properly renormalized. The D vectors 



\Lj)j' appearing in (2.1) are orthogonal by construction, and are normalized 
on their space of definition (the left part of the system, i.e. sites to the left of 
i); these shall be the only relevant vectors in the left part of the system which 
will contribute to the full analytical expression of the global state. Now, the 
trick of the trade, is considering that was obtained from the reduced 
density matrix of the previous DMRG step pf_^ = Qj which 

of course was renormalized qj = 1. This means that, the D states \Lj)f_-^ 
joint together with the local levels \s) at site i, are enough to generate the 
set of \Lj)j^: 

D d 

k=l s=l 

Here Ajf^'^^ represents the decomposition over the product basis; it can be 
either understood as a three-indices tensor (indices being s, /c, and j), or. 
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since typically D ^ d, a cZ-long array (through s) of square D x D matrices 
(from j to k). Preservation of orthonormality among \Lj)^ states determines 
a condition upon A; indeed, assuming that the local basis |s) is orthonormal 
by definition, one finds that the transformation must satisfy the equation 



D 



ki ^kj 



■'1,31 



(2.3) 



k=l 



where the superscript * stands for complex conjugation. Eq. (2.3) can be 



rewritten in an even clearer form once the A are intended as D by D matrices: 

d 

1, (2.4) 



with ■ being the standard rows-by-columns matrix product. Equation (2.3) 
follows directly from the fact that 



D 



5. 



k,m s,t 



(2.5) 



but {s\t)e = 6s,t by assumption, and {Lk\Lm)\ 



5k.m is the inductive 



hypothesis, thus (2.3) 



Moreover, looking at (2.4) under a quantum information perspective, we 
clearly understand that the ^If' actually form a set of Kraus operators for a 
completely positive trace preserving (CPT) map [1]; CPT maps are the most 
generic transformations mapping density matrices into density matrices, they 
represent the action of a quantum channel on an open system (for details, 
see appendix |A). On our case, the set of ^i^' define exactly the CPT map 
A^CPT performing the inverse DMRG transformation — )■ p^_-^ as follows 



(2.6) 



s=l 



Let us now go back at (2.2); as DMRG procedure is recursive, one can 



apply the same argument several times, for instance, until reaching the first 
site. This leads to 



1^.)^= E (^S -4? ■•••■4!') 1^1)1® 1^2)2®... (2.7) 

S\...Sl = l 

the component of \Lj)^ over the local homogeneous product basis made of 
\si. . .sg) is now expressed in terms of a product of matrices. These matrices 
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have a number of rows and columns always bounded by D, although in 
general it is impossible to require for all of them to he D x D square matrices 



and satisfy (2.3) at the same time; typically, as the left boundary grows near, 



their size shrink, up to the first site, whose A^s} are all one-row matrices. 
Similarly, Ai^^ are one-column matrices. 

In fact, we can associate a correlation space dimension Dfj to every site 
to the left of i (more appropriately: to every bond i'), and state that the 



matrices AsJ have common size D£/ x D^'+i. In order for (2.3) to hold, the 



following inequality is a necessary condition: 

De <d- De-i {2.i 



for i' < i, where it is intended that Dq = 1. As (2.8) provides an upper 
bound to the correlation dimension, so does the parametric renormalization 
dimension D, typically acting as a simple cutoff. Indeed, in standard DMRG 
algorithms it is a natural choice to adopt Di/ = mm{d^ , D}. 



2.1.1 Completing the picture: 
single center site DMRG 

So far we understood how to represent in a clear, simple analytical way of 
representing the left block states \Lj)f of our DMRG. In order to complete 
the picture to include the whole system we first need to identify which specific 
architecture of DMRG (of those proposed in literature) is being used. For 
simplicity we start with the case where the DMRG optimization is performed 
by considering at every step a single center site (system) and the left and 
right blocks (environment), as in ref. [1]. 

According to such description, the D left environment renormalized states 
joint with the right environment renormalized states \Rj)i', and the 
site levels \s)£, generate the DMRG state of the whole system: 

D d 
j,k=l s=l 

the components tensor C^^]^ defines uniquely the precise state within the 
DMRG space. 

Now, the very argument we used previously to prove that for \Lj)f_^ 



eq. (2.7) holds, can be applied in a similar fashion to right environment 



vectors as well. Precisely, if pf = Pj \Rj)f{Rj\ is a density matrix 
obtained by recursive renormalizations starting from the right boundary (site 
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L), we have that 



R 



E ■ • • • ■ • ^£0 1^^+1)^+1 ® • • • ® 1^^)^- (2-10) 



Si+i...SL = l 



Like before, we encounter a product of matrices BgJ (£' > i); but notice that 
this time the preservation of orthonormahty property for the \Rj)f states goes 
from the right boundary towards the center site, i.e. propagating toward the 

[£'1 

left. This means that the matrices BiJ should satisfy a relation which is 



different from (2.4), namely: 



1. 



(2.11) 



s=l 



To satisfy the present equation the following constraint on matrices dimen- 
sions {Bi ' being a -D^'_i x Di' complex matrix) is due: 



De<d-De+i, V£'>£; Dl = 1. 



Consistently, the CPT mapping associated to Bg performs the inverse RG 
transformation, i.e. towards the right 



(2.12) 



(2.13) 



After all these considerations, we can put (2.7) and (2.10) into (2.9) 



the state I^E'dmrg) appears automatically expanded in the natural separable 
basis: 



where Cf} has been written as an array of matrices as well. Notice that the 
term within the parentheses is a scalar due to the fact that the /^l}} matrices 
are actually row vectors (one-row matrix) and the are column vectors 



one-column matrix). Equation (2.14) tells us that all the components of 



|\E') over the canonical basis are the product of L matrices, which are local 
objects and depend only on the state of the site they are associated with. 
This is the definition of Matrix Product State El. 



Once we require that matrices ^i^^ and Bt^ respectively satisfy (|2.4|) and 



[£'] 



(2.11) in order to preserve orthonormahty of environment states, the proper 
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normalization of the state becomes an equation involving the element 
Cl^ alone. Indeed we can explicitly calculate the norm of |^) by exploiting 



its MPS representation (2.14), as 

d 



Sl...SL = l 

which, after some algebraic manipulation, reads 

Si.. .3^ = 1 

. . . ^[^-iiEfiEtm^ti^-ii . . . fiti^+Hcti^Ufii . . . AtfUtW] . (2.16) 

sl-1 Sl — J- Sf+i s^ S£„i S2 si \ ' 



And, by exploiting (2.4), (2.11) and cyclicity of the trace, all the As and 



Bs matrices disappear from the equation. In the end, we are left with 



D 



1 = {^\^) = J^Tr [clfCtg = Y: (2-17) 

se=l j,k=l s=l 

the desired normalization condition. We would like to remark that manipu- 



lations performed in order to derive (2.16) are identically suitable if we were 



to calculate the one-site reduced density matrix of |\E') at site i, namely 



a 

,J = Tr,c [|vl/)|(vl/|]= J^TrbMcfj \s){t\ 



s,t=l 



(2.18) 



where, clearly, the partial trace spans the complementary of i. Simi- 
larly, the reduced density matrices pf and pf are easily accessible one the 



representation (2.14) is at our disposal. In fact, those read 



D 



D 



E E^'i^*Si^')(^i' (2-19) 



j,k,'m=l s=l 



j,k,m=l s=l 



and all the other reduced density matrices achieved through the original 
DMRG algorithm can be generated starting from the previous expressions 



via (2.6) and (2.13). In practice, the MPS representation (2.14) provides us 



a quick access to the whole information of the DMRG algorithm, and at the 
same time is more immediate and flexible than DMRG itself, proving a useful 
computational tool as we will see later on. 



2. 1 . MPS CONSTRUCTION FROM DMRG 
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2.1.2 Double center site DMRG 

The original DMRG protocol proposed by White [2] , and most of the DMRG 
architectures still in use nowadays adopt a slightly different picture than the 



one we presented in (2.9). The basic idea is to consider the active system 
block on which to perform the minimization as it were composed by two 
adjacent sites in stead of just one, coupling to the left and right environments 
as before. Of course, at fixed renormalization dimension parameter D, this 
procedure is more expensive from a computational point of view, but provides 
a big gain in algorithm precision and rapid convergence; moreover, it allows to 
manipulate symmetries in a more natural and flexible fashion, thus improving 
algorithm stability. 

D d 

I^dmrg) =Y.H ^ii ® 1^)^ ® 1^)^+1 ® \Rk)f+i, (2.20) 

j,k=l s,t=l 



In order to recover a complete analytical expression of the form (2.14), some 
manipulation on the components tensor Tj^* has to be made. The simplest 
path to take, is to consider two composite indexes a and (3: a representing 
the pair {j, s}, while /3 representing This allows us to write Tap as a 

matrix from index a to /3, of dimension dD x dD. At this point we perform 
a Singular Value Decomposition (SVD) upon T: 

Tap = A\^^ \, B^;'\ (2.21) 

where A^^^ and ijl^+^l are unitary matrices, and the diagonal matrix of singular 

values is positive semidefinite, i.e. > V7. If we write again A'^l and 

in the original j, s and fc, t indices it is clear that they satisfy the proper 



orthonormalization propagation requirements, respectively (2.4) and (2.11). 



This tells us that (2.20) can also be interpreted as follows 

dD 

I^dmrg) = J] ® |i?^)f , (2.22) 

7=1 



where, following the formalism of (2.2) we substituted 

D d 



(2.23) 

D d ^ ^ 



k=l t=l 
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Indeed, as \L^)^ (resp. \R^)f) form a set of orthonormal vectors for the 
left (right) partition of the system, equation (2.22) actually represents the 
Schmidt decomposition of I^I/dmrg)? cut at site i. The Schmidt coefficients 
must satisfy the normalization condition 1 = (^'l^') = J2-y -^7! they are the 
positive square roots of the probabilities p in the (dD-renormalized) reduced 
density matrices of either partition of the system. The latter read 



dD 



7.L 



7=1 



dD 



(2.24) 



7=1 



In conclusion, the SVD decomposition (2.21) allows us to recover a matrix 



product expression, substantially identical to (2.14). Precisely 



1*)= E (-^S 



. . AM . ijf +11 ..... ijm 



\si...sl), (2.25) 



si...sl=1 



where A'^' is intended as the diagonal matrix with elements X'iyK Normally, in 
order to press further on with DMRG algorithm, the left (or right) density 
matrix should be properly renormalized to be D x D dimensioned; but this 
is straightforward, by just cutting off the smallest singular values A.^ until 
only the D largest of them remain, and renormalize as follows 



A. 



D 



ie{l..D} 



(2.26) 



so that state normalization (^'I^E') = 1 is preserved. Even after this cutoff, 
Als] will still satisfy the condition ( |2.4[ ): this descends automatically from 
the fact that if any number of columns are cut out of a unitary matrix, a 
left-isometric rectangular matrix {A'^ A = 1, but AA^ = P = = P"'' 7^ 1) is 



satisfy (2.11). 



obtained, thus (|2.4|). Similarly, the B^^l^ resulting from the cutoff will still 



It is trivial to make (2.25) formally match (2.14), we can either identify 

— yllf' A'^1, or alternatively Cl^'*'^' = X^^^B^^'^^ and recover the previous 



a 

formalism. Similarly, we can manipulate (|2.14l) to appear in the latter form. 



such operation will be clearer once we introduced a state-invariant transfor- 
mation (gauge) of the MPS representation, which we are going to review in 
section 12.51 



2.2. VALENCE BOND PICTURE 
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2.2 Valence bond picture and 
MPS entanglement 

Interestingly enough, it is possible to interpret Matrix Product States in a 
way [6l [To] that clarifies many of their quantum correlation properties, often 
referred to as valence bond picture. The basic idea is to start by considering 
an auxiliary system space that is actually larger than the proper Hilbert 
n = c'*^ of our ID (open boundary, so far) system; then we project onto the 
original system by means of a local transformation. Let us associate to any 
site of our quantum chain a pair of D-dimensional spins, one per bond formed 
by that site, the starting state is prepared so that every pair of virtual spins 
corresponding to the same bond, is initially in a maximally entangled state 
|$+) = \aa), known in literature as entangled bond. Then apply a 

local on-site map 

^''^ = EE4i»^0'^ir (2.27) 

s=l j,k=l 

to every site i G {1..L}, where is a canonical state in the local phys- 
ical space at site i while \j, k)f^^ is a vector of the (double spin) respec- 



tive auxiliary space. Equation (2.27) applied on the initial valence bond 



state (8)^^[^1((8)£' leads to an expression where auxiliary indexes 

of neighboring A^j^^ are contracted. Then, by writing any tensor A^j^^ as a 
set of d complex D x D matrices, the state we are describing is naturally 
expressed in the matrix product form 

l^)= E (AW.Ag.....Ait:i-4?)ki---^L)- (2.28) 

Sl...Si = l 

In general, not only the A'^^'^ operations, but even the auxiliary dimension D 
of the entangled pair |$+) can be site dependent; this way the A^^'^ matrices 
are Dg^i x Di dimensioned (where Dq = Dl = 1 to ensure that the complete 
Matrix Product expression is a scalar quantity). 

It is important to focus on the fact that, since the A^^^ are basically a 
LOCC transformation (i.e. achievable by means of Local Operations and 
Classical Communication) its action can only degrade entanglement, thus 
the entanglement of the resulting state |\E') is bound by that of the initial 
state, which is known and straightforward to calculate. Precisely, consider 
the entanglement entropy related to a left-right partition of the state |\E'), 
say at bond {£,£ + 1}. This is by definition the Von Neumann entropy of 
the reduced density matrix to the left (or right) part of the system, and it is 
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bounded by the entanglement of the original pair across the bond: 

5vN (pf) = -Tr [p[ log p^] < log De; (2.29) 

where log D is the entanglement of a maximally entangled spin pair of dimen- 
sion D, like |$+) (to check this, just consider that p£ = Tr[|$+)($+|] = j^t, 
thus 5vN(pg) =\ogD). 

In conclusion, a Matrix Product State, i.e. quantum state on a ID lattice 



allowing an analytic representation as in eq. (2.28), has well-defined upper 
bounds on its entanglement. The entropy related to a left-right partition of 
the system is bounded by the logarithm of Di, with Df being the dimension 
of the Matrix Product bondlink i we are breaking. 



2.3 Completeness of Matrix Product 
State representation 

The previous observation involving entanglement in MPS becomes even more 
meaningful once we will provide a theorem of completeness of MPS represen- 
tations. Indeed, we are going to prove that, as long as we are NOT imposing 
a finite bound to the maximal MPS bondlink dimension D, any ID finite 
lattice state can be expressed exactly as an MPS. 

The argument behind this claim is quite simple indeed. Let us choose a 
site i within the (open boundary) lattice, 1 < i < L. Let l^') be the global 
quantum state, and let us consider the Schmidt decomposition of |\E') where 
the first subsystem is made by sites {1..^} and the second by {i + 1..L}: 

l^) = i;AL'l|^«)[®|i?«)f- (2.30) 

a 

Following the formalism of previous sections, iLa,)^ are left block Schmidt 
vectors and \Ra)f the right block ones. Aa are the Schmidt coefficients 

(Sa Aa = 1), but now the number D^ of values the index a can assume is 
not anymore defined a priori; instead it depends of the specifics of the |\E'), 
precisely on its partition entanglement across the bond {£,£ + !}. Similarly, 
we could adopt the same argument when partitioning the system between 
sites i — 1 and i, namely 

De-i 

\^) = Y.^^t'^\^-)li®\R.)f-v (2.31) 

a 
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Now, since both descriptions are exact, and the fact that the block {1..^} is 
actually the composition of block — 1} with the site i alone, we must 
conclude that the set of product states of the form \La)f_^i ® generate 
every \La)i state (completeness argument). In fact, we may define the de- 
composition tensor A'^l as follows 



A 



a,l3 



L 



(2.32) 



so that we can expand \La)i in equation (2.30) in the new product basis 

1^) = E E E K'l^ 4") ® 1^)^ ® 1^")^ (2.33) 



a=l 



=1 s=l 



Of course, the completeness argument we used poses a relevant constraint 
upon dimensions of Schmidt decompositions; in particular as \La)f are or- 



thogonal, they are linearly independent, and since the \La)i_ 



\s)i basis 



can generate them, it must be that De < d ■ -D^-i. Then, by construction, if 
we write A^^l as a set of matrices (from /3 to a) then it holds 



(2.34) 



s=l 



where the positive diagonal matrices A^^^ are given by a|^' 



^a,is - SaA^'f?^ and 
correspond to the Schmidt-basis reduced density matrices of the partition, 
i.e. A'^I = Pi = pf- The previous equations resume together the orthonor- 
malization preservation relation (2.4), and the CPT mapping propagation of 
reduced density matrices (2.6). 

Similarly to (2.32), one can perform the formal expansion into site ^ and 
reduced environment for the right block of the partition, where we can define 



B 



R 

/3k 



\Ra) 



R 



(2.35) 



which allows us to write, provided the completeness constraint upon Schmidt 
dimensions -D^-i < d ■ De holds. 



Di-i Dp 



1^) = E E E (^^'' ^S) ® 1^)^ ® i^")^""' (2-36) 



a=l /3=1 s=l 
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and of course, complete positivity relations read 

^ (2.37) 
^5tM.A[^-i].5M=AM. 

In the end, by applying recursively either the left-block or right-block argu- 
ment presented in this section, we are allowed to build the analytical MPS 
representation of the original state. 

In fact, for any given state and any choice of £ {1 < i < L), one can 
formally express it as 

l^)= E (4^.... AM. A[^K5gt;l.... -5^)1.1....,), (2.38) 



Sl...SL = l 



Where the matrices A^J and BsJ are respectively given by (2.32) and (2.35); 
they are -D^'_i x Di^ dimensioned, and are well defined since the Schmidt 
decomposition exists for any partition of the system. Since the choice for 
the site i to start from, is completely arbitrary, the constraint on Schmidt 
dimensions holds both left-ways and right-ways for every site in the lattice, 
namely 

De-i < dDp and < dD^-i Vf, < i' < L, (2.39) 

where, of course, Dq = Dl = 1. This concludes the proof. 

It is now important to point out a major fact concerning the completeness 



of MPS representation; since the dimension constraint (2.39) are quite weak, 
if the state we are dealing with has no limitations on its entanglement prop- 
erties (which is the typical case for, say, a random state in the many-body 



Hilbert) such MPS representation is poorly efficient. Indeed, (2.39) tells 
us that the largest correlation dimensions ii'max = max{De} are typically 
reached next to the middle of the ID chain: precisely we have 

De<mm{d^,d^-^} — > D^^^<d^/'^. (2.40) 

Therefore, in general, the typical dimension (number of rows and columns of 
Ais}) of the MPS representation does scale with the full size L of the system, 
and in the worst case scenario it grows exponentially. 

This is the main reason why, in literature, when speaking of Matrix Prod- 
uct States most of the time one actually refers to the manifold of quantum 
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states allowing an MPS representation for which the maximal bondlink di- 
mension D is finite, does not scale with the system size L, and is typically 
small. By putting together the valence bond picture (introduced in section 



2.2) and the completeness argument, one can conclude that a Matrix Product 



State representation of bondlink D can describe exactly any state whose par- 
tition entanglement is bound by logD. Equivalently, every finitely- correlated 
state is a MPS. 



2.4 Area Law and successfulness of ID MPS 

After such preliminary considerations, interpreting matrix product states 
as variational tools turns straightforward. The typical problem we want 
to address is finding the ground states of a given, typically short-ranged, 
Hamiltonian upon an OBC system with L sites. Similarly to the DMRG 
procedure, we choose arbitrarily a maximal bondlink dimension D allowed for 
the simulation that should lead us to the ground state itself, and regard the 
elements A^s^ in MPS representation as variational tensors/matrices. Then, 
we adjust variational parameters according to some algorithm (see section 



2.9) in order to minimize the energy. 

Due to the completeness theorem of MPS representations, we know that, 
for any global system size L it exists a finite D for which the exact ground 
state is representable by a D-bondlinked MPS, and such D is related to the 
estimated entanglement e of the state itself, like e ~ logD. But now we can 
exploit some theoretical knowledge involving ground states of many-body 
systems, known in literature as the area-law of entanglement [TJ El E] : The 
partition entanglement in a ground state of a non-critical local Hamiltonian 
scales with the surface of the partition itself, and not with the parted volume. 
For ID non-critical systems, this means that e does not scale with the size of 
the system L, but rather saturates to a finite value. This also suggests that 
the bond dimension D required to achieve good precisions in representing the 
ground state does not scale with L. In practice, for many tested models, the 
D necessary to get an outstanding approximation to the GS is surprisingly 
small, regardless to system size [iA\ . This very argument allows us to address 
even problems with a large number of sites, and yet deal with them in a quasi- 
exact fashion. Of course this also explains the great success of DMRG for 
ID non-critical systems, a reason which was not yet fully understood in the 
'90s. 

Indeed, the area law argument also suggests that finite-D MPS should 
also be capable to characterize a ID problem directly in the thermodynamical 
limit as e converges to a finite value (we will discuss this approach in section 
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3.4). A special interest within this framework is raised by critical ID systems 
|33j . They are known for violating the area law of entanglement by a loga- 
rithmic (with L) correction to the partition entropy, with the proportionality 
constant given by the central charge C of the model |i331 135] : 

e{^) = 5vn(pi..l/2) ~ ^ log L + C" (2.41) 

and therefore the appropriate D to represent the ground state faithfully, does 
scale in the end with the system size, according a power-law like behavior 
D oc L*-^/^ where the exponent is C/6 itself. Now since the large majority 
of the famous ID critical models have typically small central charges (e.g. 
crit. Ising, crit. XXZ, Heisenberg, have C < 1), even though D scales with 
L, the scaling function is so concave that even in that case we can address 
efficiently quite large system sizes with good precision. 

Nevertheless, it is important to remember that for critical ID systems, 
their efficient MPS representability depends directly on the central charge, 
while for non-critical systems it is natural, an automatic consequence of the 
area law of entanglement. In chapter |5] we will introduce families of varia- 
tional states more suitable to address criticality than mere MPS 



2.5 Gauge group of 

Matrix Product State representation 

By now, it should be clear that, given a quantum state on an OBC chain, its 
exact Matrix Product State representation is in general not unique. The issue 
is simple: the state components expanded in the canonical basis are compos- 
ite products of matrices, and the same product can be matrix-factorized in 
many ways. We will now define and explain the usage of a group of trans- 
formations that manipulate the set of matrices in the representation, under 
which the physical state is invariant: by definition this is the gauge group of 
MPS representation. 

Let us start again from the state whose MPS representation has 
bondlink dimension D, and is given by 

l^)= E (2-42) 

Sl...S£=l 

where matrices A^fj are -D^_i x Di dimensioned {Di < D, V£). For every 
i < L, we now define an invertible square matrix X^, of dimension Di x D^. 
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The expression within parentheses in eq. (2.42) is left invariant by 



si S2 sl 



= ■ X,) (X,- . A^S -X.)... [X,\ . Ait:' ■ Xl-.) {X,\ . Afl) . 

(2.43) 

But now, any term of the form {X^\ ■ Als"^ ■ Xi) is again a -D^-i x matrix, 
and as above we have d of them per site, indexed by s. In conclusion, the 



latter expression in eq. (2.43) is again a Matrix Product, where the bondlink 
dimensions are preserved site-by-site, and the original Matrices of the 
representation underwent the (gauge) transformation 

Af = Xi\ ■ AM ■ X, Ws e {l..d}, (2.44) 

while the state is left invariant, i.e. 

I^)= E {Bl}}.BlS.....Blt_l'.Bi^)\s,...s,). (2.45) 

Sl...Si = l 

For any nearest-neighboring bond, Xg defines an allowed transformation as 
long as its inverse is defined. Therefore, the gauge group of Matrix Product 
States is equivalent to the direct sum of the groups of Isomorphisms of Di 
dimensioned complex vector spaces 

L-l 

6;mps = 0Iso(C^^). (2.46) 

e=i 

In order to define properly ^mps we did not need to summon the Hilbert 
structure: the invertibility condition is a rank dimension requirement, not 
a metric constraint. This remark is definitely sensible, since the correlation 
spaces are fictitious, virtual, and therefore there is no reason for a gauge 
group to mingle with the physical-space metric properties. 

Finally, notice that the gauge group we built ^mps is identified by the 
initial choice of site-dependent bondlink dimensions Di, which we required 



to be left unaltered from the transformation. Actually, in (2.43) we could 
have used any rectangular x D' matrix Xi (with D' < D^) which is right- 
invertible, i.e. 

Xi ■ = 1, but . = P = P V 1, (2.47) 



and adopt such X^ in (2.44); this, of course, leaves the matrix product invari- 
ant, but the bondlinks of the representations are altered, their dimensions 
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increased {Bi now are De-i x D', and Bs are D' x D^+i). However, the 
state contains the same amount of entanglement as before, but we are spend- 
ing more resources to describe it: we are working in a non-optimal numerical 
framework. Moreover, this extension Q' to the previously defined ^mps is 
clearly a group lacking an inverse-element property. For these reasons, in 
most cases it is interesting to limit the study of MPS gauge features on ^mps 
itself, under which the MPS representation space, given by the {Di}e, is 
stable. 



2.6 The canonical form 

The presence of a gauge group for MPS provides an computational advan- 
tage, since freedom and manipulability of our description tools are increased. 
At the same time, the capability of quickly recognizing state properties, or 
comparison between states is reduced, as even MPS representations of two 
identical states may look very different, when their gauges are incompatible. 
The simplest way to avoid such difficulty is to break the gauge invariance 
by hand, i.e. by characterizing a representative in the class of equivalence 
for MPS, which is easy to achieve, recognize, and completely general. This 
concept realizes in the definition of a canonical form for MPS representations. 
We say that a Matrix Product State, of bond dimension D 

\^)= E {Am.Afi.....A['-]\.A^^)\s,...s,), (2.48) 

Sl...S£=l 

where Ag^ are De^i x dimensioned matrices, with open boundary condi- 
tions {Do = Dl = 1), is in the (right-) canonical form if it holds 

d 

J2 A^I^ ■ = 1 Wi,l<i<L 

s=l 

J2 A^f ■ A'^-'i ■ AM = AM \fe,i<e<L ^ ' ^ 

s=l 

yY[o] = j^m — g^j-^^ every A^ is positive and full rank. 

We recall that the first equation is the CPT-condition to preserve orthonor- 
mality among Schmidt vectors, when propagating from the right: for this 



reason, we will, from now on, refer to this gauge- invariance breaking, (2.11) 



(2.37), as right gauge, for brevity. 
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Given any MPS, it is always possible to write a canonical matrix prod- 
uct representation for the same state; the bondlink dimensions are equal or 
smaller than the original ones. An operational proof of this statement is 
explained in detail in ref. [10], we are now going to sketch the fundamentals, 
as many of the involved manipulations will be useful later on, and this is the 
perfect context to introduce them. 

2.6.1 Proof of canonical form generality 

Let us take an OBC-MPS representation given by Bf^ matrices 

l^)= E {Bl^}-Bf].....Blt_l^.Bi^)\s,...s,), (2.50) 

S1...SL = 1 

we are going to define explicitly a set of rectangular matrices Yi and Zg, with 
YiZf = 1, such that by applying 

= r^-i^MZ^ for 1< £ < L, 



the resulting matrices Alf' satisfy (2.49), and they represent again |\E') faith- 



fully via (2.48). Moreover, the resulting bondlink dimensions will be equal 
or sni aller than those of representation. Precisely, the full rank condition 
(2.49 ,c) for A^^l tells us that the resulting A^f^ representation uses the minimal 
bondlink dimension, for every bond, necessary to describe l^'). 

The Ye and Zi represent a gauge transformation followed by a correlation 
space truncation; constructing them is quite simple. We start from the right 
of the ID chain, by performing a Singular Value Decomposition (SVD) of 
the matrix -B]^'^ read as if Sl were an incoming index, and j an outcoming 
one: 

<=E'']>?<1.. (2-52) 

where U^^^ and A^^'^ are respectively left and right isometric matrices, i.e. 
WU = 1 and AA'^ = 1, and the diagonal matrix A'-^l is positive. Not only, 
but we can make A^^^ strictly positive, by just cutting /3 values for which 
aI,^^ = out of the sum. If we do, f/'^' and A^^^ continue to be isometries, as 
any subset of columns of U^^^ (resp. of rows of A^^^) is still an orthonormal 
set. Then we just fix Yl-i and Z^^i matrices as 

y^.i = A^^l^V^^ = [/[^IaM. (2.53) 
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By construction Yl^iB^J"^ = A^J"\ which is an isometry and thus in the right 
gauge. Similarly we define C^"^ = bI^ ^^Zi_i, and of course 

l^)= E (2.54) 

Sl...S£=l 

is still a faithful representation of the original state |\E'): the rightmost 
bondlink Dl-i might be decreased after the transformation, but due to the 
SVD argument we know we disregarded only zero components. In other 
words we could say that Zl_iYl-i = P = P^ = P^ is the projector over the 
actual support of the bondlink space (in pf^ representation). 

Now we proceed recursively: we consider the composite matrix cj^)^ (with 
£ starting from L — 1 and moving left) whose incoming index a is the pair of 
indices {k, Si}, and again we perform a SVD 

<t. = j:u5,^'^L (2-55) 

By construction A'^^ satisfies the right gauge condition, since 

E^'X"= = «^.. ^ E4^-4'S = i- (2-56) 

a Si 

Then, bondlink £ space is truncated to the support of A, and the (pseudo-) 
gauge transformation given by 

= Al^rVM^ = t/^AM, (2.57) 

which lead to = Y^.^B^^Z^ = A[^\ and redefine = CF^^ 

so that we can apply the procedure again on site £ — 1. Once we arrive at 
the left-end of the chain, ( 
initial state normalization 



the left-end of the chain, cli' is already in the right gauge by assumption of 



in conclusion A^s} = Cl^' = Bf^Ze = A^J^K This proves the first statement 
of (2.49): converting a complete MPS representation in order to be fully in 
the right gauge is operatively possible by a recursive application of Singular 
Value Decompositions. 

The second statement of (2.49) is a direct consequence of the CPT map- 
ping argument we presented in previous sections, in particular it corresponds 
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to equation (2.13). Finally, the full rank condition follows from the fact that 
after every SVD steps we truncated the bondlink space (dimension Di) to 
the support of the corresponding reduced density matrix A^^l This issue is 
argumented in details in jlOj . 

As a concluding remark to the present section, we would like to point 
out that the canonical form we just presented is the right-directed one, i.e. 
is made so that every MPS block is in the right gauge (and the correlation 
space used is minimal). Of course, we could similarly define a left-canonical 



form, where MPS is completely in the left gauge, i.e. (2.4) (2.34), and other 
statements still hold: 



J2 ■ = 1 W,l<i<L 



s=l 
d 



M (2-59) 



s=l 

m = AW = 1, and every A is positive and full rank 



such is the left-canonical form for MPS representations. The demonstra- 
tion adopted in this section to achieve the canonical form is operational in 
the sense that is exactly the algorithm we apply in numerical settings: the 
computational advantage of using canonical MPS, apart from immediate es- 
timation of entanglement, will be clear as soon as we explain how to achieve 
expectation values onto an MPS state. 



2.7 MPS and Observables 

By now, we understood that MPS are outstanding candidates as tools for 
simulating condensed matter one-dimensional many-body systems. Then, it 
is fundamental that we realize how to achieve expectation values of observ- 
ables, in a clear and efficient way. If we are to adopt, say, the global energy 
as a simulation benchmark, so that our goal becomes achieving the absolute 
energy minimum, we need first to calculate the expectation value (^|if|\l') 
of the Hamiltonian H: and operator which is nonlocal, but it is explicitly 
written as a sum of local (separable) terms. For simplicity we can assume it 
couples only nearest neighboring sites 

e=l q 1=2 p 
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Let us start from getting the expectation value over a MPS of a separable 
observable O = 0^ G'^l, where every operator ©1^1 can depend on the site £ 
on which it acts. We have 

si...s„ ri...rn 

X (Afl • . . . • (ri . . . (g) |si . . . (2.61) 

e 

Now we define the so-called transfer matrices, as follows 

e|]^ X^(r|X|.)(AM®Af), (2.62) 

s,r=l 

where X is a one-site operator, acting on site £; transfer matrices e|| are 
D^-i X -D| dimensioned. We now calculate the Eq^ for every £, and the 
expectation value becomes simply a multiplication of the whole string of 
transfer matrices 

(*|//|*) = eW - Eg] .....E^^i. (2.63) 

The computational cost for acquiring this expectation value scales only lin- 
early with the system size (recall that now we treat L and D as independent 

parameters of the MPS variational ansatz); unfortunately, there is still a 
harsh dependence on the chosen bondlink dimension D. Indeed multiplying 
two X matrices costs ~ elementary operations, yet by adopting 
some technical tricks we can further improve this scaling law: 

• We should start performing the multiplication from the right (or left) 
boundary, as \Ql-i) = Eq' is a one-column matrix, i.e. a column 
vector. And multiplying a dimensioned vector \Q,) for a x 
matrix has an overall D'^ cost. 

• instead of multiplying directly Eq |Q^) = \Qe-i) one can first calculate 
iri^l) = [Alfl ® l]\Qe), followed by |ni^^) = Es(^|0^^^|s)|ri^') and finally 
\Qi-i) = ® A'^M]|ni.^^). These operations requires respectively a 
number of elementary operations equal to: dD^, (fD^, and dD^. 

In the end. the total computational cost to achieve the MPS-expectation 
value of a separable observable O scales, with size L and bond dimension D, 
as 

#cost ~ L {2dD^ + d^D"^) . (2.64) 
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Usually, the first term in the parentheses is the leading one (and the other 
is negligible), since the typical bond dimensions chosen in simulations are 
sensibly larger than local space dimensions D ^ d. 

The result we got holds in a quite general scenario (provided that O is acts 
locally); we will now see that if the involved operator has a small support, we 
can considerably improve this limit by exploiting the gauge group of MPS. 



2.7.1 Local support Observables 



Let us assume that the observable O we are interested with does not involve 
all the sites within the ID chain, but only a small connected subset of those, 
say lattice sites between ii and £2 {1 < ii < £2 < L) . Recalling the previous 
argument involving transfer matrices, i.e. eq. (2.63) we can write 



[^1-1] 



E, 



^2 + 1] 



■ E 



(2.65) 



Notice that on the sites outside {ii..i2} the observable O acts trivially, so 
we are considering the transfer matrix of the local identity operator 1 there. 



Now, as the expectation value in (2.65) is a physical quantity, i.e. it de 



pends on the properties of the quantum state and not on its specific MPS 
representation: it is a quantity invariant under the action of the MPS gauge 
group. At the same time the transfer matrices are not gauge invariant, so it 



is advisable to choose a gauge that reduces the computational cost of (2.65) 
Precisely, we choose a gauge that turns our MPS to look like this 



.Sr=l 



. ^[^2 + 1] . 
St2 + 1 



.Bi^^)\s^...SL), (2.66) 



where MPS tensors to the left of the {ii..i2} support are in the left gauge 
C^gAlAs = 1), those to the right of the support are in the right gauge 
C^gBgBl = 1), and those in the middle can be in any gauge chosen by 
the user, with the only constraint that they must satisfy the global state 



normalization condition. As before, the form (2.66 ) can be achieved by means 



of recursive Singular Value Decompositions that define appropriate gauge 
transformations (from site 1 to ii to fix the A, from site L to £2 to fix the 



B), exactly like we did in the section 2.6 



Now we focus on the transfer matrices of the outside zone El . Consider 
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site L, due to the assumptions we made, it holds 

\Ql-i) = = E f E ^^^r BfA \jk) = |$+) (2.67) 

j,k \ s,r / 

where |$+) is the (unnormahzed) maximally entangled canonical vector 

l^-") ^ E l^'^') = E ^^M- (2.68) 



It is east to see that eq. (2.67) holds recursively for all \Qe), i > £2, since 



j,k \ s,r a,j3 / 

^5,,,.|jA;) = |$+), (2.69) 

where we used both the fact that the operator acts like identity on site i 
(thus Sg^r), and the recursive hypothesis \Qi) = | $'*'). An identical argument 
can be applied to transfer matrices to the left of the support of O, where left 
gauge condition can be exploited to see that {Qi\ = ($^| for any i < ii. The 
conclusion simply follows: 

(^|0|vl>) = ($+|Eg;|....-E£||<|.+), (2.70) 

which means that the number (T^cost) of elementary operations we have to 
perform does not even scale with the full size of the system L, but merely 
with the size of the operator support: 

#cost ~ (£2 - ii) {2dD'^ + dPD"^) . (2.71) 

Honestly, we traded the modest effort of performing the SVD, needed to 
convert the MPS in the proper gauge, to obtain a faster (and non-scaling) 
computational speed in acquiring finite-range MPS physics. 

We will see that this result can be partially exploited even when we are to 
compute expectation values of observables which are not local and not even 
separable, but allow a natural decomposition into local terms, such are the 
Hamiltonians of typical short-range interacting models. 
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2.7.2 Hamiltonian-like Observables 

We are now interested defining an operational algorithm that, exploiting 
MPS properties, computes efficiently the expectation values (\E'|iJ|\E') of an 
operator H which is formally written as a nearest-neighboring Hamiltonian 



of the system, i.e. like (2.60). For algebraic reasons which shall be clear soon. 



we rewrite it as if = H^, where 

L 



= E E 0? + E E H'' ® 0"? ' (2-72) 

£=£' + 1 q 1=11+2 p 

As before, it is important that we focus on the computational cost of this 
data acquisition. We learned that working in the proper MPS gauge is in- 
strumental for economy of calculus, thus we already start from a canonical 
MPS representation (say the right one) 

l^)= E W5g...5ff)ki...^L), (2.73) 

Sl...Si = l 

where all the are right-gauged, requirement which also guarantees proper 
state normalization (^E'l^E') = 1. 

Again, our scheme to acquire (\E'|if|\E') has a recursive formulation: we 
need to propagate the contraction of our MPS structure and at the same 
time include every term of H. Since this is expensive by definition, we will 
try to regroup and sum the partially contracted tensors every time we can. 
Then the starting point, at the right boundary, is defined as follows: 

j,k \ q r,s 

\xt'') = E (thPirlQ'fl^) ] Bf^B^f^^ljk). 

j,k \ r,s 



(2.74) 



Now we propagate towards the left. The idea is that l^^^^) shall contain all 
the elements of the Hamiltonian who have support in {i + 1..L}, while |xp^) 
will take care of nonlocal neighboring terms across the bond {i — Of 
course, we can exploit the gauge conditions, telling us that \Qi) = |$^). This 
allows us to calculate every \xp^) directly 

IX^') = EE (th?(r\Q'f\s)] (2.75) 

j,k a \ r,s / 
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Instead j^'^ is obtained via the recursive relation 

j,k a \ q r,s / 

+ EE (j:t.('-\^'>)] Bj?BtSb*)("/3|xl?) + 

j,k a,0 \ p r,s / 

+ EEE«;b^)Mi^l')- (2-76) 

j,k a,0 s 

By means of the transfer matrices formaUsm, we can rewrite the two previous 
equations in a more compact, and clearer, form 

1 p (2.77) 

Acquiring all these data, for every £ requires an overall computational cost 
(apart subleading trends) of 

#cost ^ L{i^q + 2#p + 1) {2dD^ + (fD^) , (2.78) 

where #g and #p are respectively the number of one body and two body 
terms in the Hamiltonian expression. As you see, even in this complex sce- 
nario, the scaling behavior with D and L of the computational cost remains 
roughly the same. Now we can conclude that 

(*|^ri*) = EW-ESf'-...-Ef]|e[^l), (2.79) 

and in particular |\E') = \^^^^) which is a scalar number, and it is exactly 
the energy of the state if H is the actual Hamiltonian of the system. 

The great improvement in computing expectation values we encountered 
for separable observable is recovered in case of (short-range) Hamiltonian 
operators: the full computational cost to acquire the energy (which will be 
later adopted as variational functional) scales only linearly with the system 
size. If one thinks that (full-search) exact methods typically bear an expo- 
nential cost in L, it is easy to understand why DMRG/MPS architectures 
are regarded with great interest. 
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2.8 Pictorial representation of Matrix 
Product States Tensor Network 



Through the present chapter, we learned how to deal with MPS: their math- 
ematical properties that allow algebraic manipulation (gauge group), and 
their physical properties that allow to control quantum entanglement. On 
the other hand, the equations we encounter start to look cumbersome and 



confusing, like eq. (2.76). To work around this issue, we are now going to 
provide an alternative way to express MPS representation that is based on 
diagrams and graph theory rather than standard analytical expressions. This 
will prove a faster and clearer fashion to represent states, observations, and 
matrix multiplication; and will become instrumental in later chapters. 

Let us start back from our definition of MPS. If our system is a ID-OBC 
lattice with L sites, and {|s)}s is the local canonical basis, then a generic 
state of the system is written as |\E') = Tsi...sl\si ■ ■ ■ ^l)] the complex 

tensor (with L indices) Tsi...sl uniquely defines l^'). Now, stating that |\&) is 
an MPS (with fixed bondlink D), is equivalent to say that Ts^...sl allows the 
following decomposition 



De<D 



■11 /i[2]s2 /i[3]s3 
31,32 32,33 



.A 



[L-1]SL-1 a\L\sl 
'3l-2,3l-i 3L-1 ' 



(2.80) 



As mentioned before, all the A^ elements are three-indices tensors (apart the 
first A'^I and last A^^'^ MPS blocks, which have only two indices due to open 
boundaries), being the physical index, i.e. related to the local canonical 
state Is)^, while the two being the correlation space indices linking to the 
two neighboring MPS blocks, namely A^^"^' and A'^"*"^!. Equation (2.80) tells 



us that a MPS is the result of a multiple contraction of (possibly variational) 
tensors, or more simply a Tensor Network. 



Let us draw eq. (2.80) in the following pictorial graph 



AW 




Aid— 


A[3] 






1 




1 1 



_ • • • A^^M 



where every block (graph vertex) represent a single tensor A^^\ the legs/links 
attached to it being the indices s (vertical one), j^.i and (horizontal 
ones). Connecting two tensors though a given link means contracting the 
product of the two over that index. By these very simple rules, one sees that 
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(2.80) is recovered, but there is no need anymore to write down either the 
sums or every single index s or j by its name. Everything is imphcit in the 
pictogram. 



As we discussed in section 2.5 , on each of the contracted index/connected 
hnks we can insert an isomorphism Xi together with its inverse X^^. Since 
their global action cancels out during link contraction, they have absolutely 
no effect on the global tensor T. This is exactly the gauge of Matrix Product 
States, which we represent like this 







^"^/^ I : 




1 






1 1 




AM 


^^^^P'' ■ 1 




m ' 








I 1 




1 


r 1 



meaning that B\ 



[e]s 



,.a^a,k and B.,^ 



[i+l]s 



s 

a,k' 



This oper- 



ation is definitely equivalent to (2.44), once we have defined an isomorphism 
X^ for every link i, 1 < i < L. 

To break by hand the freedom granted the gauge group, we defined two 
particular gauge choices, namely the left and right gauges. We said that 



AM is in the left gauge if T^s^s^s = 1, i.e. 
equivalently, read in terms of transfer matrices. 



E 



s,a j,a ■'^a,k 



IE. 



pictorial representation, it becomes an equation between graphs: 



or 



and in 





(2.81) 

Similarly, the right gauge condition is the left-right specular of this graph 
equation. Having such equation in graph form lets us to see immediately 
how to exploit the gauge by substituting pieces of the Tensor Network and 
thus eliminating tensors. In particular, assume we are to calculate the state 
square norm (^|^) = E{.,} "^i-^^'^i-^l = 
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If all the MPS tensors are in the left gauge, we can substitute (2.81) into 



the diagram, and tensors start to literally cancel out, from left to right, until 
only 1 = (^E'l^E') remains. Moreover, the left gauge condition ensures (2.6) 
(reduced density matrix propagation via CPT map), which corresponds to 




AW , 






A 





(2.82) 



In section 2.7.1 we saw that having the leftmost MPS tensors in the left 
gauge and the rightmost in the right gauge provides a great advantage when 
computing observables having local support. Precisely, let O = (R)!?-*, B^^l, 



and we require that ^^^4''"^ Alf' = 1 for £ < 
i > £2- Then, the part of the graph outside { 



and ^^Ai'^Utfl = 1 for 
£2} cancels out thanks to 



(2.81) and we are left with 
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= ($"*'|Eq^' ■ . . . ■ Eq^' I'^'"'^); as we saw in the previous section. 

To explain the algorithm we use to calculate the expectation value of 
a nearest-neighbor interacting Hamiltonian H, we introduced the relations 
(2.76) and (2.77), where the transfer vector was defined via a recursive 



scheme. Despite the unclear look of those equations, it is possible to resume 
them both in a simple and intuitive graphical equation: 





where we just regrouped together the one-site operators as R[ 



(2.83) 



and the two-site ones: R^^^ = J2p ® © p • 

Several equations that we will encounter in this thesis involve Tensor 
Network contraction or decompositions, in most cases they allow a diagram- 
matic version, granting immediateness, and clarity of understanding. So, 
where appropriate, we shall provide it for completeness and comfort for the 
reader. 



2.9 Minimization algorithms 

We mentioned MPS being powerful variational tools for simulating ground 
state of ID many-body systems, and we also managed to give a prescription 
for evaluating the energy of the matrix product state. It is finally time we 
adopt such energy (\E'|if|\E') as a functional for variational simulation, and 
describe an algorithm that drives our trial Matrix Product State toward the 
absolute minimum of this functional. Despite the huge reduction in varia- 
tional parameters we are left thanks to the MPS representation (~ LdD^ 
rather than d^), performing full search of the minimum within the whole 
parameter space at once is still too expensive for practical purposes. Instead, 
we will follow a scheme similar to original DMRG: the idea is to perform only 
local or quasi-local variations of the whole state representation, namely vari- 
ating a limited number of connected MPS blocks while keeping the other 
fixed. Then we repeat, while choosing each time a different compact subset 
of MPS blocks to variate, until convergence is eventually reached. Like tra- 
ditional DMRG algorithms, usually one or two adjacent block are variated 
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at a time, and refrain by sweeping towards the left or the right (bouncing off 
the boundaries of the system, if our problem is OBC). 



2.9.1 Single variational site 

In this framework, at every minimization step only one block of the MPS is 
being treated as variational, say the one related to site i; the other ones are 
fixed, and in practice the energy functional itself will depend on them. We 
also assume that all MPS blocks to the left of i are in the left gauge, and 
those on the right are in the right gauge, so that 



l^)= E {A^!-----^l--^-C!i-Bltl'-...-Bi^)\s,...s,), (2.84) 

Sl...Si = l 

with At^'l (E.^^r^i''^ = 1) and (EsBf'^B^V = 1) Axed, and we 
are searching the C^^^ which minimizes {'^\H\"^). It is always possible to 
gauge transform an MPS to achieve form (2.84) by means of repeated SVD 
as we saw previously, the whole singular part of the decompositions has been 
embedded inside C'^l 

It is easy to calculate the explicit dependence on C^^l of the energy: by 
means of transfer matrices formalism one can write 



+ (2-85) 



Let us interpret it: the first and second term embed respectively the terms of 
the Hamiltonian having support to the right and left of i; the third term are 
those acting on i only, and the two last terms couple i with its neighbors. As 
all the terms contain just one E^^l, the functional is quadratic in the tensor 
CM: = E. C*m,n'^mnrClfc> whcre the effective Hamiltoman H is 

hermitian. Also, we should take into account state normalization = 
E. ^*],kCj,k^ which is a constraint of the problem, and therefore must be 
inserted in the functional with its appropriate Lagrange multiplier e. 

In conclusion, the resulting Lagrangian reads £{C,C*) = (\E'|iJ|\E') — 
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= {{C\n ~ eAf\C)) 



c = 




s 




(2.86) 

where |C)) is intended as a dDf^Dn^i ~ dimensional vector, % and A/" as 
observables on such space, representing respectively the effective Hamiltonian 



and effective square norm. Here "H is given by the same five terms of (2.85) 
in the same order of appearance they read 




(2.87) 

with 1^^), as before, obtained recursively via (2.83), and similarly (^^~^| arises 



from the left-right specular equation. 

You see that, thanks to the gauge condition chosen for the A^^^ and 
B^^'\ the effective square norm M coincides with the identity operator on 
the (iD^-dimensioned effective space, so that £(C, C*) = {{ClTi — £'^\C)) as 



prescribed by (2.86). Finding the minimum of a quadratic Lagrangian is now 



straightforward since, by exploiting the fact that C and C* are independent 
for complex differential calculus, one has 



dC{C, 



|0)) 



'H\C))=e\C)). 



(2.88) 



Therefore we have to deal with a standard eigenvalue problem for Ti, and 
among solutions 'H\C^ = b\C)) we have to consider the one giving minimal 
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value of ((C|?i|C)) = e{{C\C)) = e; in other words, we have to find the lowest 



eigenvalue solution of the problem (2.88), and we are done. In conclusion, we 



mapped an eigenproblem for the whole dimensioned global system into a 
dD^ eigenproblem for the local tensor C'^'. When Z) is a parameter chosen 



by the user and not dependent on L, solving the local problem (2.88) requires 
the same effort for every system size. Also, the great advantage of working in 
the proper gauge framework is clear now, since if A/" were not to coincide with 
1, we would have to deal with a generalized eigenproblem ("HlC)) = eAf\C))): 
more expensive and less stable. 



It is clear that after one has found the minimal |Cmin)) for (2.88) the 
energy of the resulting |\l/2) is necessary decreased (or equal) from the initial 
guess 

(*2|i?|^2) = ((C^minl^lC^min)) < ((Cg^ess I^H | Cg.ess)) = (^ll^l^l)- (2.89) 

This is how the algorithm proceeds towards energy minimization, after this 
step we move to the MPS block to the immediate right i+1 (or the immediate 
left i — I) and repeat. Of course, to complete the iteration step one has to 
perform the proper gauge transformation so that, after finding C^^^, turns it 



into left gauge — )■ A^^l (or right gauge if we are sweeping left) so that (2.84) 



immediately holds for site i + But this is easy, just perform a singular 
value decomposition of C^^]^ 

cSk = J2<P^pU^^>^ (2-90) 

with a being the composite index {j, s}, and A'^A = UW = 1. The gauge 
transformation is then 



(2.91) 



which concludes the iteration step. 

The algorithm is usually carried on until some convergence threshold in 
the energy e has been achieved. In most cases this simulation procedure 
converges surprisingly fast, as very few sweeps are necessary to reach a sta- 
ble minimum, even if we were to start from a completely random variational 
MPS. In several simulations where the single-site framework was adopted, 
computational results are in good agreement with theory and/or experiment, 
yet this protocol presents some difficulties. The monotonicity of the energy 
functional at every iteration step, even if it allows fast convergence, hides 
the possibility of getting stuck in local minima of the variational parameters 
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landscape: in order to work around this issue one has to insert manually 
artificial fluctuations, as proposed by S. White in DMRG context Sim- 
ilarly, the algorithm encounters trouble when dealing with symmetries (see 
appendix |B]) , where the user is forced to insert symmetry-breaking fluctua- 
tions by hand. 

Despite on how we can solve, with big or small success, these issues in 
the single-block framework, a very common way to work around them is 
to recover the original idea that gave birth to DMRG, i.e. dealing with a 
two-site block minimization at once. 



2.9.2 Double variational site 

This time we want to variate two adjacent blocks at the same time, say 
C^^l and C'^"'"^', while keeping flxed the other ones. The most clever way to 
do this is forgetting that C'^' and C^^~^^^ are two distinct MPS blocks: we 
consider them as a single overall tensor Mj^,!^'''^^ = Cf^^^ C'f fc*^^ on which 
the Lagrangian functional is quadratic, and adopt M as our only variational 
element. Of course, this allows us to momentarily describe more entangle- 
ment across the bond {£, £ -t- 1} than what would be normally allowed by a 
D-bondlink MPS. Therefore, to provide an iterative scheme, we will embed 
in the algorithm a method for entanglement truncation, so that in the end 
the original MPS representation is recovered. 

Then, let us start from our initial guess for the iteration step 

Sl...SL = l 

(2.92) 

where every M^^^s^+i is given by the matrix product cf} ■ Cl^+i^ and MPS 
blocks As J (resp. BsJ) are in the left (right) gauge. As a whole, M is 
a tensor with four indices, and notice that alongside the sudden increase 
of allowed entanglement (from logD to log dD), an increase of variational 
parameters, with respect to the standard MPS case, comes out: from 2dD^ 
to d^D\ 

As we mentioned the Lagrangian is quadratic in M, and thanks to gauge 
relations for A and B matrices, the effective normalization M is again the 
identity operator, since (^'l^') = Tljksr l^/fc P- Then 



C{M,M*) = {{M\n\M)) -e{{M\M)), 



(2.93) 
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n 




+ 



+ 



c 




(2.94) 

Graphs 1 and 2 contain the terms of the Hamihonian having support outside 
+ 1}; graphs 3,4 and 5 those of support inside {i,i + 1}, and the last 
two represent interaction of the inner sites with the environment. 

Like for the single site case, the optimal M is found via lowest eigenvalue 
problem, i.e. smallest e allowing 



n\M))=e\M)). 



(2.95) 



Now we have to manipulate the newly found M in order to recover the 
standard MPS form. 

To do this we proceed again via singular value decomposition. First we 
write the tensor M as a matrix Ma/3 where the composite index a ~ {j,se} 
refer to the left bondlink index j and the left physical index sg, while /3 ~ 
{k, S£+i} to the right bondlink k and site S£+i indices. Now we calculate the 
SVD as 



dD 



M, 



a/S 



7 ^7/3 



dD 

7 



-^7 -°7fe 



(2.96) 



so we can stick it into (2.92) and obtain again the MPS representation. Still, 
the bondlink i has increased dimension D"*^™ = mm{dD£_i, dDg^i} ~ dD. 
But now the positive values are the Schmidt coefficients of the partition at 
bond i of the state l^f), and (^|^) = XI 



dD 

7 



)i2 



Therefore, in order to recover 
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the best approximation of this state allowing \ogD entanglement, we must 
truncate the smallest Schmidt coefficients A^, until only D of them remain. 
Namely, if were sorted in decreasing order (A^ > A^+i), we keep only the 
first D of them, and renormalize to preserve state norm 



A. 



A, = ^=^ ie{l..D} -> E^? = l; (2-97) 



moreover, A^ and 5^ will satisfy respectively the left and right gauge con- 
dition even after the truncation. 

We can now write, assuming we are sweeping towards right, M^^ ^^^^ = 
A^i ■ CT,:? , where C^i+^l^^^^ = ~X,B^^:t'^'^^' so that 

S\...SL=\ ^ ' 

(2.98) 

The MPS representation is now ready to perform the next algorithm itera- 
tion, just identify the new two-site block Mg^^^^s^^^^ = Cl^+i^ ■ -Si^+a'- 

The double-site based algorithm we just described presents two impor- 
tant improvements with respect to the single-site one. First, as two adjacent 
blocks are being modified at the same time, reconstructing the correct short- 
range physics runs much faster, and since the Hamiltonian is made of nearest 
neighboring terms the energy is extremely sensitive to the n-n physics (es- 
pecially for non-critical systems) thus leading to a faster minimization con- 



vergence. Secondly, at the time we perform the truncation (2.97), we allow 
for errors in our description, as we force the state to carry no more entan- 
glement than the MPS representation allows. Therefore, slight fluctuations 
appear, identified by eventual small increases in energy. This is actually an 
advantage of the protocol, as fluctuations are a natural way to discourage 
the algorithm from getting stuck in local minima of the energy landscape. 



2.10 Matrix Product Operators 

So far, we applied the Matrix Product formalism in order to build multi- 
indexed tensors 7^i...<i^, which were components of a target state |\E') over a 
separable (canonical) vector basis \si . . . sl)] but it is clear that its capabil- 
ities extend to every algebraic construct which can be expanded in a basis 
of local elements, regardless their nature. No doubt, applying the Matrix 
Product concept to describe nonlocal operators seems the most natural goal 
to pursue. 
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Matrix Product Operators (MPO) were firstly introduced in [36] and 
used, for instance, to describe eitlier tliermal mixed states, time evolution 
paradigms for MPS [M] [SH], or long range interaction Hamiltonians 
|10] . Their open boundary formulation is definitely similar to that of MPS: 

d d 

= E E (^ff" ■ • • • ■ ^ff •• • • • • ^^1' (2.99) 

SI...SL ri...rL 

where, if \se) is the canonical vector basis on site i, then \se){re\ is the canon- 
ical local operator basis. The local-basis expansion is then performed over 
these canonical elements, while adopting Matrix Product-based coefficients. 
To every site i, incoming physical index r^, and outcoming physical index sg 
we associated a Dg^i x Di matrix Of}'^\ which sums altogether to a four- 
indices tensor on every site: 



L JL JL _l 



om 



r ^ 



0[3] 

T 



(2.100) 

As you can see, blocks 1 and L of the MPO have only one correlation space 
link index to be consistent with the OBC setting. As for MPS, Matrix 
Product Operators are typically prescribed according to a maximal bondlink 
dimension D [Di < D V£, with D non-scaling with L) which makes the 



expression (2.99) manageable for practical purposes even for large system 



sizes L. Such D also poses a limit in the entangling capabilities of G, actually 
binding the amount of long-range correlation the operator can create. 

MPOs are outstanding tools when the goal is to apply a transformation 
O to a state |\E') whose MPS representation is available. In fact, the resulting 
state 6|\E') is automatically expressed in Matrix Product form: 



0i^)= E (or -...■OS-) (4;' ■...■4 

{se}{ri}{te} 

x\si...SL){ri...rL\ti...ti) = E (4^----^l?) ki---^^^) (2-101) 
where 



(2.102) 
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Am 

























Bin 



1^ J ■ 



Truly, the bondlink dimension of the target MPS G|\E') is increased to D' = 
Da ■ Do, the product of the original MPS bond Da and that of the MPO 
Do- So, it looks that application of MPO to MPS is definitely expensive in 
terms of the bondlink. This is true, and nevertheless easy to work around: 
it is sufficient to reduce the target MPS to the desired D" properly. This is 
quickly done by following the usual steps: 

• Choose a bond, say + 1}. 

• Gauge transform the MPS so that blocks to the left (resp. right) of 
the chosen bond are in left (right) gauge, so that Schmidt coefficients 
of the partition emerge explicitly 

• Truncate the smallest Schmidt coefficients and renormalize to one the 



remaining ones (squared), as (2.97) 



• choose another bond and repeat, until every bondlink has been renor- 
malized to D" or less. 

The error we intake when renormalizing the state is compatible with the 
amount of entanglement we are discarding (which is explicitly known by 
comparing Von-Neumann entropies before and after truncation). Being able 
to transform MPS into MPS becomes fundamental, for instance, if we want to 
describe a time-evolution of a system whose starting point is a finitely corre- 
lated state: within this paradigm it is very useful to understand how to write 
an MPO representation of a given Hamiltonian, and how to exponentiate it 
efficiently. This is a major point of interest of ref. 



2.10.1 Matrix Product Density Operators 

A relevant class of operators we are typically interested in, is the family 
of density matrices, i.e. positive, unity trace, operators. Although it is 
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instructive and useful decomposing such operators into MPO form, it is even 
more interesting to exploit their positivity (as well as positivity of any partial 
trace, i.e. degree of freedom reduction) to further decompose their Matrix 
Product structure. Indeed, if we consider that any p > can be written 
as p = XX"^ (and, conversely, the whole space of matrices X generate the 
class of positive operators via XX'') one is encouraged to build the MPO 
decomposition of X rather than p itself. Doing so not only eliminates the 
positivity restraint on the resulting MPO, but also gives us an edge for dealing 
with state transformations, as the application TpT"!" simplifies into TX, which 
could be a nontrivial numerical improvement. 

Precisely, in ref. [36] the Matrix Product Density Operators (MPDO) are 
properly defined. They are those MPO, according to ( 2.99[ ), whose blocks 
are given by 



(2.103) 



r=l 



where di is at most dD^^iD^. Decomposition (2.103) is actually splitting the 



matrix product layer into two sub-layers stacked together, as 





(2.104) 

where the pentagonal shape of the tensors in the diagram specifies that ten- 
sors in the upper layer are up-down specular to those in the lower layer (plus 
complex conjugation). 

If we were to give an interpretation to the tensorial index dimension 
we could invoke again the valence bond picture: indeed log di is the maximal 
allowed entanglement that the system can share with an external degree of 
freedom coupling expressly to site £, e.g. a local thermal bath. 

MPDO are useful tools for addressing one-dimensional open systems, es- 
pecially where the mixing with external media acts on the bulk itself. It is 
even possible to formulate master equation problems with matrix product 
formalism. 
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2.11 Example: exact Matrix Product State 
representation for Slater Determinants 



We would like to conclude this chapter with a simple, yet practical ex- 
ample of the Matrix Product formalism applied analytically to a specific 
class of many body states: Slater Determinants. In fermionic problems, 
Slater Determinants are the starting point of most many-body calculations 
(like Hartree-Fock); they are states where N fermions share no quantum 
correlations, each one of them filling an orbital which is typically solution 
of the mean-field Hamiltonian. Nevertheless, since such orbitals are non- 
necessarily localized in a chosen configuration space, they can still manifest 
self-correlation entanglement in the separable basis. So it is probably the 
simplest among non-trivial matrix product decomposition problems. The 
following construction is somehow related to ref. |10], but I developed it dur- 
ing my Philosophiae Doctorateship as an independent project, supported by 
G. Santoro and V. Giovannetti. 

Let us deal with spinless fermions, for simplicity: the first step we have to 
perform is to match the physics of this context with the algebraic formulation 
adopted so far in this chapter. To do this, we can completely forget about 
physical dimensionality of the problem and boundary conditions: the only 
initial structure we need is a complete set of L one-body wavefunctions. We 
also choose a complete ordering for these. They will represent the sites in 
our ID OBC (spin) system, placed according to the chosen ordering; the 
canonical local basis corresponding to |0) = empty level and |1) = filled 
level. The mapping of a fermionic system into a spin system is made via 
standard Wigner transformation 



where \Q) is the vacuum state, q the destruction operator on level i, and a 
being Pauli matrices. Now any state |\E') can be expanded in such product 
basis l^') = J2s-^ sl '^si...sl\si ■ ■ ■ sl), and we are going to apply the matrix 
product formalism to the components tensor T^^...^^, local dimension d = 2. 
For sake of completeness, let us even write the explicit MPS expansion in the 
second quantization formalism. 



ivi/)p_= j2 {A^!;-...-Ai'i){c\r...{c[r-\n). (2.106) 



where the construction operators cj are placed in the correct order, and 
obviously (cj)° = 1. Simple as that. 





(2.105) 




2 



Sl...SL = l 
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A Slater determinant state |E) is defined as follows 

|s) = £l4...4|fi) 



(2.107) 



where cj^ fills a one-body orbital which may have a nontrivial expansion over 
the original one-body levels we chose as basis. Precisely the transformation 
is given by 



(2.108) 



0*(£) being the first-quantization decomposition of the Slater orbital a onto 
the original wavefunction basis i. Orthogonality is required among a orbitals, 
i.e. ^£ 0o(^)0^(^) = Sa,i3. As we know from literature, such considerations 
lead us to write |S) in its explicit determinant form 



0<£i<...<£]v<i 



01 (4) 



01 (^2) 
02(^2) 



01 (^JV) 

02 (^iv) 



0iv(^l) 0Ar(^2 



Ml 



cj ...c{ Jfi). (2.109) 



To find a MPS representation for will be instrumental to give a 
Matrix Product Operator description for Fermi operators Cq, over delocalized 
orbitals. Since the vacuum is already in (trivial) MPS form |0 . . . 0) we will 
then find the MPS structure of |S) by stacking together MPOs of cj^ according 



to (2.107), as we did in (2.101) and (2.102) 



2.11.1 MPO for delocalized Fermi operators 

We are now going to provide a MPO representation of cj, which is compact, 
elegant, and very general. The only ingredient we need is knowing the ex- 
pansion of the orbital a in the original one-body wavefunctions basis 0a (^). 
Then the goal is finding the i?'^' satisfying 

1 1 




E E (&o|fii^-...-5if^|6L)ki....L)(ri...r^ 



si...Si=0 ri...ri=0 
1 

= E (&o|i?i^■...■sif^l&L)(c^)-...(4)^-l^])(^]|c2-...c-, 

{s},{r}=0 

(2.110) 

where we explicitly set vector boundaries to the matrix product expression 
(actually |6i) is a vector and (60 1 a functional) so that we will able to define 
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every Bg homogeneously, even those at the furthest sites. In particular we 
need D = 2, and the solution we found is given by 



1 
1 



1 



bK 














aii) 



0aW a- 



1 
-1 



|0) 



(6o| = (0 1)=(1| 



(2.111) 

where the information on 0q,(£) is used on only one of the 16 elements of the 
four-indices tensor B^^'^ 
homogeneous in 



Apart from that the expression (2.111) is formally 
as we wanted. 



To show that (2.111) reproduces the correct action of cj^ we first notice 



that when performing the matrix product contraction, the terms which con- 
tain one and only one a~ are the ones that survive: indeed (1|0) = (l|cr^|0) = 
0, while (l|(j~|0) = 1, but on the other hand (cr~)^ = 0. So we can reduce 
the expression (2.110) in a simple sum over the site £ upon which the a 



IS 



activated, becoming 



a 



£-1 



(2.112) 



which proves the equivalence. As an additional remark, it is easy to see that 
it is possible to deactivate the global action of such MPO by just changing 
a correlation boundary state, say (6o|- Precisely, if we were to set (6o| = (0| 
instead of (1|, the MPO expression (2.110) would coincide with the identity 



1 instead of cj^; one can then regard the correlation space boundaries as local 
switches that control the whole matrix product behavior, even if it is not 
localized. 



2.11.2 MPO stack to MPS 



We can now go back to the Slater determinant state |S) = c\cl . . . cjy 
and adopt the engineering we learned to define its whole MPS exact repre- 
sentation (2.106). In particular we start from the vacuum state |00...0), 
which is trivially an MPS with D = 0, and we apply cjy to obtain again an 



MPS. Then we refrain, by applying in order cjy.i, c 



and so forth, up 

to c[: every step is performed following the prescription of eq. ( |2.102[ ) (al- 
though we should never renormalize if we want our description to be exact). 
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In conclusion we have 



(ji...(jjv=0 



5 [^,2] 93 



92 



91 



(2.113) 



As you see, our description uses as a whole a total bondlink dimension of 
D = 2^, regardless from L. Actually, since every Bq is the null operator, we 
can also restrict the previous sum to qi > q2 > ■ ■ ■ > Qn, since every term 
for which any < qk+t, with t > 0, would give zero contribution. In the end 
it is a sum of merely N terms. Similarly, we define the correlation boundary 
vectors: 



\bL) = |0) 



( 1 \ 





and (6o| 



(0 ■■■ 1). (2.114) 



Putting these ingredients together leads us to the decomposition of our Slater 
Determinant S in the MPS representation, where explicit boundaries of the 
matrix product expression are present 



(2.115) 



Sl...SL=l 



Let us briefly analyze the matrices we built via (2.113). It is easy to see that 



Aq' it is always the identity = Idxd, while J^^ contains the information 
upon orbitals, expanded in the original wavef unctions. To make this clear, 
we show as an example the cases = 2, for which it holds = 0i(£)[l ® 
(y~\ ^ (\)2[t)\(y' ® a% i.e. 



N 



1 











o\ 




























\ 





-02 (^) 


01 (^) 


0/ 



(2.116) 



It is clear that, the only matrix products that lead to nonzero amplitude are 
those where two excitations |1) are present. Thus the sum ( 2.115[ ) reduces to 



0l(^l)02(^2)-0l(^2)02(4 



(2.117) 



where we have recovered explicitly the determinant expression. Also, let us 
write down the case with three orbitals to be filled A^ = 3, in this scenario 
we and up with 
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/ 00 00\ 
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3(r) 
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Mr) 





-Mr) 


Mr) 






the reader is invited to check that the resulting amphtudes are correct. 



2.11.3 Efficiency of the description 

We want now argument that, if no further information upon the orbitals (pai^) 
being filled is exploited, the exact representation we just gave is the most 
efficient in terms of MPS. By this statement we mean that we are spending 
the smallest bondlink dimension needed to faithfully reproduce the correct 



amount of correlation the state can manifest. From section 2.2, we know 



that a D-dimensioned bondlink MPS allows up to iS < logg D entanglement, 
i.e. Von Neumann entropy of a partition (the logarithm base of 2 is chosen 
as common ground in quantum information theory), thus the optimal is 
D = 2'^. Now, inequality S < N is guaranteed by the existence of an exact 
MPS representation (2.115). But if equality S = N can be achieved for some 
choice of (f>a{^), we also proved representation optimality. 

To obtain it, we just adopt a special set of (doubly-periodic) plane-waves 
0a (^) = explAniai / L) . For simplicity let us perform a half-system partition, 
and define a new double set of orbitals {<f)\t\i),(j)a^{i)}a from the previous 
ones as 

0[^](£) = v^e(L/2-£)0„(£) and 
<j)l^\i) = V2Qii-L/2)Mi), 

with being the Heaviside step function. Even though in a general case a 
new set of wavefunctions generated via (2.119) would no longer be orthonor- 
mal, it is clear that with the specific choice of L/2 periodic plane waves, 
orthonormality is preserved: ^i(j)a'\^)(f)*^^\i) = 0, as the supports are dis- 
joint, and ^e(j)a\i)(l)*f\i) = T.i<P^<^\^)(P''f\^) = ^a,i5- So we can define 
Fermi operators corresponding to this new set, satisfying the anticommuta- 
tion rules {ca,„Cfs^,} = {c„,l, c],,/?} = and {ca,L,cJ, i} = cj, j^} = Sa,(s- 

It is clear that the original c decompose in the new ones as Ca = 2^^/^(ca,L + 
Ca^n), thus letting us write the whole Slater determinant state as: 

1^) = i {^i,L + • • • (4,L + 4,k) \^), (2.120) 
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Of this state, we want to calculate the density matrix reduced to half the 
system, say the left one, so we trace over the right-half degrees of freedom 
Tr^ [|S)(E|]. With this goal, we set = Ci|S) and consider: 



Tr 



R 



R) 



2 v<LTr«[|S')(Sl]5i,L + Tr«[|S')(S'|] 
2l2x2®Tr^,[|S')(S'|] = ^®pf, 



(2.121) 



where we used the cyclicity of the trace over right support operators Ca,R, 
and clearly ci^r\T,') = 0; then we noticed that p^' and {c\j^pY ci^l) have 
orthogonal supports. Now we repeat the same argument on |S'), and proceed 
by induction. In conclusion, we can claim that is (isometrically equivalent 
to) 2-^1 2iVx2Jv, the maximally mixed state on a 2 dimensioned space, whose 
Von Neumann entropy is just N. This concludes the proof. 

An intuitive, but not naive, interpretation of such result can be given 
in the following terms. As fermions occupying the various orbitals must be 
mutually uncorrected due to the Slater determinant state nature, the only 
possible entanglement the system can manifest under a real-space partition 
is given by the self-correlation of every orbital, separately accounted. In 
fact, in the studied case, we presented uncorrelated completely delocalized 
orbitals, each one carrying the entanglement of a unit (i.e. the amount of 
entanglement shared by a spin singlet), so A^ is naturally the total amount. 



OVERALL REMARKS 



The (T^ matrix in equation (2.111) is the one and only responsible for 



establishing the correct anticommutation relations of Fermi statistics. 



That said, it is straightforward to modify (2.111) so that the corre- 



sponding MPO is describing a Bose operator instead: you just need 
to replace B^^-^ = 1 and leave the rest unchanged (also extensible to 
abelian anyons by using phase gates B^^\ = e^'-'^"'). 

We mentioned that the present design is modeled on spinless fermions, 
but actually is naturally extensible to fermion with spins. The only dif- 
ference is that at the very beginning, when we are selecting a complete 
basis of orbitals, we need to specify a complete basis of spin-orbitals 
instead, and then choose a complete ordering. Any ordering is fine 
and does not compromise the MPS cost in terms of D as long as the 
particles are uncorrelated and fixed in number. 
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2.11.4 Tensor grid representation of one-body 
wavefunction basis change 

Let us recall that, when we derived the MPO representation for cj^, we also 
mentioned that it is possible to control its overall action by adjusting the 
left correlation boundary vector (6o|: namely the MPO coincides with cj^ if 
(6o| = while it is just the identity 1 for {bi\ = (0|. In other words 
1 

J2 ■ • • • ■ Bl'^^^-^\b,) cr . . . cr\^){^\ cl' ■■■c?= ci\ (2.122) 

{s},{r}=0 

with q G {0, 1}. Also recall that the fermionic orbitals 0q,(£) we filled to build 
the Slater determinant state were an orthonormal set: let us complete it to 
an orthonormal basis {0a(^)}a, with a G {1..L}. The dimension must be 
L by the assumption that the original set of L wavefunctions was complete. 



For any of those 0o(£) the corresponding MPO is given by (2.111). 

Now we stack together the MPOs, like we did for the Slater state, but 
instead using only of them, we stack the complete set, ordered from a = 1 
on top to a = L at the bottom; moreover, instead of using the standard left 
correlation boundary vector (6oU = (1| we set a generic (6o|a = {lal- It is 
obvious that the operator arising from this construction is equivalent to 

gt9i gt52 gt93 ___gt^te_ (2.123) 

Finally, we apply such operator to the vacuum \VL). The meaning of all this 
construction is that we actually defined an application on the binary strings 
of {qa\a to the real Fermi space, as 

(gi...gz.|-^5l'^5^.--£^l^)- (2-124) 

By linearity, this map extends to all the space generated by ((?i • • • Q'l I , which 
corresponds to the whole correlation bondlink space (as (q'l • • • (Jl | is its canon- 
ical product basis). The map is clearly bijective and thus invertible. But 



you notice that the inverse of (2.124) is basically a Wigner transformation 



from the Fermi space to its spin representation where this time the (pa have 
been chosen as basis of one-body wavenf unctions, so it is formally similar to 



(2.105), but the basis has changed (the old one is associated to the c, the 
new one to the c). 

In conclusion, we could use all this MPO stack formalism to represent a 
many-body state transformation |^) corresponding to a change of the chosen 
basis of one-body wavefunctions. I.e. assuming that we can expand 

{4=0 {q}=0 
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then the two components tensors T°'*^ and 7"'^'^™ satisfy the equation: 






r T T T 



(2.126) 

where the blue tensors in the grid are exactly ^I^-"] of ([2TTT1), with i being 



the coordinate in the horizontal axis, and a the one in the vertical axis (the 
origin is the lower-left corner). The upper and rightmost edge tensors are 
trivially |0) and |0). 

An interesting remark to this result, is that the grid Tensor Network 



that appear in (2.126) can be efficiently contracted, despite having several 



closed loops in its geometry (see section 4.5), thanks to particle-conservation 
symmetry relations. 



2.11.5 Extensions to Configuration Interaction 

In quantum chemistry settings, the simplest path to move beyond the mere 
mean field paradigm is adopting Configuration Interaction. In those de- 
scriptions, Hartree-Fock solutions are adopted as a canonical vector basis of 
orbitals for further calculations. According to such viewpoint, one is inter- 
ested to express correlations by superposing few to several Slater determinant 
states, which typically share some of the HF orbitals as well as differ for 
other ones. If the energy minimization problem were to be performed over 
the whole space of Slater states the result would be exact, still this would be 
an extremely hard problem: thus generally the amount of orbitals for which 
the involved Slaters differ, is kept to a small, manageable number. 

Having this scheme in mind, we would like to extend our previous Slater 
MPS (MPO stack) representation to embed also Configuration Interaction 
states, where different orbital excitations are coherently added. The ultimate 
ingredient of this perspective would be writing Matrix Product representa- 
tion for every operator generated by the Fermi ones c}^ through sums and 
multiplications. Of course the related zoology is huge, so we will limit our 
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discussion the simplest nontrivial case. 
Consider for instance the expression 



©2+2 = a c\cl + /3 clcl 



(2.127) 



We want to describe O2+2 as a Matrix Product Operator, and as you can 
guess there is no unique way to perform the extension from the normal Fermi 
operator case. Depending on whether we focus on the adaptability of the de- 
scription or the economy on the bondlink dimension we end up with different 
proposals. 

Standard Guess - this path exploits the standard technique to sum 



coherently Matrix Product objects, and is strongly based on (2.111); thus is 
highly suitable for further generalization, but at the cost of a sub-optimal 
bondlink dimension. Let us now adopt D = 8 and consider 



r[^,2+2]' 
j 















(2.128) 



where the tensors are those defined in (2.111) for cj^. The basic idea 



behind this construction is to use a correlation space which is the Cartesian 
sum of the two original correlation spaces, and a matrix product object 
which is the block diagonal composition. Similarly we define the correlation 
boundary vectors, which contain information on a and f3: 






/3 
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*V ) 




\0j 



(2.129) 



where we used distributivity of the tensor product ® with respect to the 
Cartesian sum ©. Similarly, (60 1 = (11)® (0 001). 

Cheap Guess - this path focuses on keeping the lowest correlation 
bondlink dimension possible, and actually requires D = Q. 



J-6X6 
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and 



/ 1 











o\ 





-1 

















-1 














-1 

















-1 





Vo 














(2.130) 



while boundaries are as before = |0) and (6o| = (o . oi) = (5|. By 
multiplying the i^I^'^+^l matrices it is easy to see that we are reproducing the 
correct action of the operator, i.e. 



+ 



J2 {« {m^i)M^2) - M^2)M^l)] 

+ /3 (M£i)MQ - 04(^2)03(^1)) j 4 4- (2-131) 



Like previously, we argued if this Matrix Product representation is optimal 
in terms of correlation bondlink dimension: we found that a state of the form 
©2+2!^) has a real-space partition entropy of entanglement at most equal to 
5/2. This implies that a faithful MPS description would require a D > \/32, 
so that D — 6 is the smallest allowed integer, and thus optimal. 

The present proposal presents various options for generalization, although 
finding the analytical MPO expression for a generic operator which is cheap- 
est in terms of D is definitely a hard task. With this last speculation we 
conclude this analytical example of Matrix Product formalism for interesting 
states in condensed matter physics and quantum chemistry. 



60 



CHAPTER 2. MATRIX PRODUCT STATES 



Throughout this chapter we dealt uniquely with open boundary condition 
problems, and developed a formalism of Matrix Product States based on 
the OBC framework. Of course, such a description can be adjusted to fit 
naturally periodic boundary conditions as well, taking care of the correct 
amount of entanglement. In the next chapter we will introduce a periodic 
description for finitely-correlated states and thus Matrix Product States, with 
its proper formulation and tricks of the trade; this will be instrumental in 
the proper definition of a thermodynamical limit. 



Chapter 3 



Periodic and infinite Matrix 
Product States 



One of the major issues for standard DMRG architectures in ID problems 
is deahng with Periodic Boundary Conditions (PBC). It was soon clear that 
traditional DMRG ideas could not be applied to PBC with the same success 
and simulation precision, but it was only with the advent of MPS representa- 
tions that this trouble become clear and argumented. Indeed, while in OBC 
the DMRG procedures describes all and only the finitely correlated states, 
i.e. those states whose entanglement is bounded by a finite value (which typ- 
ically does not scale with system size L) in PBC the correspondence is not 
exact any longer. Nevertheless, finitely correlated states play again a very 
important role in describing ground states of short-range interacting models, 
as they manifest the correct entanglement area-law. Indeed, even in PBC 
finitely correlated states naturally lead to a matrix product representation, 
but the formulation |4H H2] is slightly different from their OBC counterpart. 



3.1 Valence bond picture for Periodic MPS 



In section 2.2 we introduced the valence bond picture to argument and con- 
textualize MPS with open boundaries; its is straightforward to extend such 
description to a periodic system. To every site we associate a pair of spins, 
each one D dimensioned {D chosen by the user, often sensibly larger than the 
local degree of freedom dimension d). We prepare this virtual state so that 
every pair of neighboring sites share a maximally entangled state through 
the D-dimensioned spins l^"*") = D~2 \aa) (entangled bond). Notice 
the difference with the OBC case, where we had L sites and thus L — 1 phys- 
ical bonds: in PBC every site has two neighbors (there is neither first nor 
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last site, or, if you prefer, sites 1 and L are neighbors), so the amount bonds 
is L. The virtual-to-physical mapping is defined identically to the OBC case: 



D 



(3.1) 



=1 j,k=i 



As before, which we are going to apply it to the composite entangled bond 
state <^iA^^^{<^f, |$'^)^/"^/+i). Immediately, one can see that the resulting 
state can be expressed as 



(3.2) 



where the Ai^' (resp Ai^-^) are no longer vectors (dual vectors) in the corre- 
lation space, but matrices, ^ -Dlso < D, like for every other site i. The 



trace operator in (3.2) makes the inner matrix product cyclic, so there is no 



starting nor ending point of the ID ring. Also let us represent |^') 



AW 




(3.3) 

diagrammatic version of (3.2). If we now are interested in estimating the 



entanglement of a connected subset of sites, we can use the same argument 
for OBC and get a similar conclusion. In fact, if we want to part the system 
in an given interval of sites and its complementary, we need to break two 
entangled bonds of the virtual state And since the resulting 

state |\&) has entanglement bounded bounded by the first one, we have 



^VN [p. 



= -Tr[p£,...£2 logy 



< log D,, + log D,, ~ 2 log D, (3.4) 



which is twice as in the OBC case, where we could split the system while 
breaking just one entangled bond. 

An interesting point concerning periodic systems is dealing with transla- 
tional invariance symmetry. As most models have translationally invariant 
(TI) Hamiltonians H, exploiting the expected TI of the ground state be- 
comes fundamental for every simulation method. This is mostly true for 
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PBC, where TI is meaningful and spontaneously broken only in exceptional 
cases (when ground space degeneracies arise), whereas in OBC the breaking 
is naturally induced by the presence of boundaries. 

It is immediate to see that if the tensors in the MPS representation do 
not depend on the site, i.e. A'^l — y A regardless from £, then the state 
is translationally invariant: 

d 

J\^)= Tr[As,-As,-As,-...-As^]J\si...SL) = 

si...sl = 1 

d 

= Yl Tt[A,-A,^-...-A,^-A,,]\s2...slSi) = \^), (3.5) 

si...sl=1 

where T is the elementary translation operator. The original state is obtained 
again by using trace cyclicity and a relabeling of the indices s. A more press- 
ing problem is the inverse: given a translational state T|\l/) = allowing 
a periodic MPS representation, does it have also a homogeneous representa- 
tion, i.e. where matrices are not site dependent? We will constructively, and 
positively, answer such question right away. 



3.2 Translational MPS admit a 
homogeneous description 

Assume we are starting from a site-dependent MPS representation A^^ of a 
state as in (3.2), we will build another MPS rep. B where matrices do 

\ 

(3.6) 



not depend on the site any longer. Let us write 

/ 



B, = L-L 



A 



[2] 



^[^-1] 
V / 

we will now show that the MPS built with these matrices is equivalent to the 
original one. In fact 



Y Tr [B,^ ■ ■ . . . • B,j] |si . . . sl) = 



Sl...SL=l 



q=0 si...SL=l 
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Figure 3.1: Construction of the homogeneous MPS B (3.6). the resulting 
bondhnk can be seen as composed by the channel carrying j G {1..D} and 
that carrying q G {^--L} so its total dimension is LD. the starting MPS A^^l^ 
is now seen as a four-legs tensor, as the site label i becomes a tensor index. 



The black dot is the triple Kronecker delta, while A 



9,9 



5n' 



g' ,g+l mod L • 



L-1 



L-1 



E E 

q=0 si...SL=l 



Tr 



At 



■•••■41, \s^...s,) = -J2v\^) = \^), 



q=0 



(3.7) 

because T|\E') = |\E') by hypothesis. In conclusion, we succeeded in building 
a homogeneous representation for a generic finitely correlated state on a ID 
PBC ring, but not without expenses. Notice, indeed, that the bondlink 
dimension we end up with is D' = Yle-^i ~ -^-^5 with Di being the original 
bondlink dimensions of A'^l; a linear scaling law with the system size L arises. 
At the same time, we are describing the same amount of entanglement as 
before, so the B representation is definitely sub-optimal. 

Unfortunately, this is common ground when dealing with PBC Matrix 
Product descriptions (for both states and operators). Regarding this issue, 
ref. [lO] proposes a trivial example involving the W-state T'|0 . . . 01), 
which has minimal MPS bond dimension of 2, that necessarily increases to 
L if we want to give a homogeneous MPS description. 



3.3 Expectation values for periodic MPS 

Often in numerical settings, dealing with periodicity is more tricky and ex- 
pensive than corresponding OBC versions of the same problems: for MPS 
this peculiarity manifests immediately in an increase of computational costs 
needed to acquire expectation values of observables. 
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In section |2.7| we defined a procedure to evaluate physical quantities of a 



MPS in OBC with a number of operations that scales nicely, as (2.71), with 



main simulation parameters: system size L and correlation bondlink D. In 
the best case scenario, where the separable observable O acted on an interval 
of £ sites, we estimated the cost to scale like ~ £D^d. 

Scaling laws are not so nice for periodic MPS, and the compactness of 
operator supports does not he lp, d ue to the presence of a global closed loop 
in the graph (see also section 
and we are looking for the expectation value 



4.5). Precisely, assume that O = (3)^1^ 0^^' 



X 



si...Sri ri...r„ 



X Tr 



^ n • • • ^ r-i 



(ri...ri|(g)eM|si...Si). (3.8) 



We can still rewrite this equation in a simpler Matrix Product form thanks 
to the formalism of transfer matrices, defined identically as before e|| = 
^^^(r|X|s)(y4[f' ® A*^^^). The difference is that this Matrix Product is also 
cyclic, i.e. 



(^|0|^) = Tr 



IL0 ■ ILi 



■ E 



[^1-1]' 



(3.9) 



Multiplying two transfer matrices costs ~ D^, an expense that can be reduced 
to 2dD^ + (fD^ by calculating in the order = {As ® 1)E, then = 
X;^,r(r|X|s)M^, and finally E' = Xlrl^ ® ^r)Qr. Unfortunately, this is the 
only improvement that can be made in general. 



Equation (3.9) has no right and left boundary vectors, which were instru- 



mental to remove a scaling power out of the cost. Moreover, the gauge 
group can be no longer exploited to eliminate terms from the product of 
matrices; this can be argumented as follows. We would like, for instance, to 



transform the Ejf' into the identity so that it disappears from (3.9). However 



the MPS gauge group transforms the transfer matrix E^' according to 



E^: 



(3.10) 



But the input matrix Ey could be entangling, and the transformation is local 

[£1 

and invertible, so there is no chance that a generic Ey can be mapped into 
a non-entangling operator (like 1) this way. 

In conclusion, if we want to acquire the exact expectation value of a 
product observable on a PRC MPS, the computational cost is 



#cost ~ L {2dD^ + d^D^) . 



(3.11) 
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Honestly, it is absolutely convincing that for large L the system will be less 
sensitive to finite size effects, thus manifesting an emergent physics which 
is very similar to OBC physics in the bulk. We could exploit somehow this 
limit to reduce costs while acquiring controlled errors; nevertheless, it is 
useful to understand how to work with MPS in the thermodynamical limit 
before elaborating this idea. 



3.4 Thermodynamical limit MPS 

The chance of extending a Matrix Product State description so that it is 
actually representing an infinite system L — )■ oo, follows directly from the 



notion of MPS homogeneity we discussed in section 3.2 



Let Aji^ be the elementary tensor block of a PBC homogeneous Matrix 
Product State. For every system size L, this A defines a unique state in 
the 2^ dimensioned Hilbert space, thus forming a sequence of states |\E'l). 
The thermodynamical state is defined through physically relevant quanti- 
ties, namely expectation values (0)oo of compact support observables which 
should coincide with the limit of {0)l for L — ?■ oo. 



+00 



(3.12) 



Of course such expectation value limit must be well-defined in order for the 
thermodynamical state to be consistent; thus we need to understand under 
which conditions upon A this uniqueness is achieved. 

Therefore, consider the expectation value "^{O) of a compact support 
observable O, which acts on i adjacent sites 

^^O) = hm ^,{0) ^ hm \^\[ f = 



lim 

L—>oo 



Tr 


El... 


El Eo 


El 


...El 


Tr 


El... 


El 





lim 

L— >oo 



Tr 


E^' Eo 


Tr 


Ei 





(3.13) 



where Eq is the £-sites composite transfer matrix of the operator O: 



Eo = 5^ (^1 ■ ■ • ri\0\si ...s,) Tr [{A,, ® A^J ■ . . . ■ (A,, ® A*J] 



Sl ...Si = 1 

n . . .r£ = I 
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In eq. (3.13) you see that, since we are keeping A fixed, for a generic L the 
MPS state {"^l) will not be normalized, so we have to introduce manually 
the square norm in the expression = Tr[E^]. 

It is clear that the Thermodynamical limit state ^oo(") niust not depend 
on how we perform the limit itself, nor which boundary conditions we used 
at finite sizes. So we must obtain the same result even starting from an open 
boundary setting, as long as A still describes the bulk, and distance between 
the support of O and boundaries diverges. So we will introduce arbitrary 
correlation-space boundaries (6ieft| and | bright) by hand, and write 



^oo(O) = lim 



(&lcft| Eo E'^^ I bright) 



[0\eh\ li-i rrightj 



(3.14) 



we will require that this limit does not depend on (6ieft|, | bright), n or m {n 
and m being any two positive integers); it also must coincide with the limit 



in equation (3.13) 



Lemma - uniqueness of limit (3.13), (3.14) holds iff some spectral re 



quirements upon the transfer matrix of the identity operator are satisfied: 

• among eigenvalues Aq, of Ei, there is one Aq strictly bigger than all 
others in modulus, i.e. |Aq,^o| < |Ao|, 

• Aq eigenvalue is simple, meaning that only one related eigenvector |e,^) 
exists i.e. the eigenspace of Aq has dimension 1. 

Let us prove this statement. The transfer matrix E^ is not necessarily diag- 
onalizable, but as it is on complex field, we can expand it in its generalized 
eigenvector basis (which is not orthogonal in general), so it appears in the 
Jordan block form 



El 



/ Ao 1 


\ 


Ao 1 




Ao 1 




Ac 






Ai ■■• 


V 


... ) 



(3.15) 



where we highlighted the generalized eigenspace of Ao. The conditions we 
required upon Ei tells us that |Aa^o| < |Ao|, and that there is only one 
Jordan block corresponding to Aq, as it appears in (3.15). Then, given a 
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random vector \v), it is possible to demonstrate that 



with 



e^) normahzed 



hm 

q—^ca 



-^0 PO 



n\v) 



v\Ei\l\v) 



(3.16) 



1, provided that 



has no null component 
To show (3.16) first expand |f ) in the 



over the generalized eigenspace of Xq 
generalized eigenbasis |e^g), where a is the eigenvalue index, d refers to the 
Jordan block where the basis element belongs, and w its position within the 
block. Then the composite application of several eJ^ gives 



[As-v 



(3.17) 



d.w 



where V^^\q) is a polynomial function of g, of degree x, and whose coef- 
ficients depending on A; Aq is the size of Jordan block d. When we take 



the limit (3.16) of the latter expression, all the components belonging to gen. 
eigenspaces different from the Aq one vanish, as (Aa/|Ao|)'^ — )• for any a 7^ 0. 
Then, also components over generalized eigenvectors of Aq which are not the 
true unique eigenvector |e^) disappear, since their polynomial multiplier is 
of lower degree, i.e. 'P^^°~^\q)/\V]^°\q)\ — for any w 0. 



Therefore the composite action of on a generic vector \v) maps it 
(after normalization) to [cq 



that 



With a similar argument, it can be shown 
where (e^ | is the only 'left-eigenvector' (eigenfunctional) 
of Ej. Notice that by construction (e^| is not necessarily the dual of |e^) 
via Riesz representation theorem, but the two vectors are not orthogonal 
either, so that (e^|e,^) 7^ 0. After all these considerations, we can extract 
the desired conclusion, i.e. 



^oo(O) 



lim 



Tr 


E^^ Eo" 


Tr 


E^ 
''-1 





lim 

L—^ca 



E.Tr 


\v){v\Ei EoE^ 


E Tr 


|w)(«|Ef+'" 





lim 



{bleit\ Eo E™^ I bright) 



(felefti E 



{n+m)L+e 



''right J 





Eo 


e^) 






e^) 





E 





e^) 




e^) 



(3.18) 



where the limit correlation boundaries are defined by Ao|e^) = Ei|e^), and 
Ao(e^| = (e^|Ei. As you see, the result is consistent regardless if we are 
following an OBC or PBC scheme. The reverse implication in the Lemma is 
trivial by counterexample. 
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The capability of expressing the expectation value for every observable 
as in (3.18), can as well be formulated in terms of density matrices p£. That 
is, a thermodynamical limit quantum state can be properly defined by the 
sequence of reduced density matrices {p£}e for any finite size i, having the 
property that when tracing partially larger-sized ones, we recover smaller- 
sized ones: 

{pe}e ■■ — > Tri..^o [pi] = Tr^_^„..£ [pi] = pi_t^ V£o < ^- (3.19) 

Here we are implicitly considering translational invariance as well, which is 
automatically granted by homogeneity of the A in our case. Then, following 



the prescriptions of the lemma and what we learned from (3.18), we can 
write down the expression for the reduced density matrix p£ of an arbitrary 
number of sites as 



pi 



SI , 

ri , 



■ Si 





(As^^a:.^).. 






1 







\si...se){ri...ri\, (3.20) 



where the partial trace property (3.19) is an automatic consequence of the 
fact that T.r,s^rA^s^A*)\e^) = Aolco"). 

An interesting and physically relevant comparison with Matrix Product 



Density Operators, we introduced in section 2.10 can be made once we gave 



the pictorial representation of (3.20): 



Pe = 




(3.21) 

where ^ = [AQ(e^|eQ*)] ^ is put for correct state normalization, and correla- 
tion boundary vectors are the solutions of the eigenproblem 




A* 




AnX 





A* 



(3.22) 
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for maximal modulus eigenvalue |Ao|. The resemblance of (3.21) with (2.104) 



is evident, in particular our MP-thermodynamical state is written in a pe- 
culiar MPDO form where the mixing dimension di = 1 (meaning no exter- 
nal degrees of freedom that couple locally) except for the boundaries where 
di = D. Moreover, the two designs match exactly iff e^, read as a matrix 
(being a two-indices tensor, consider one of the two as incoming matrix in- 
dex, and the other one outcoming), is positive. In this case we can write 



XXK and embed the X within the MP-block A^!^ Bf = A,X as it 



was a gauge transformation, and we find precisely the formalism (2.104). 

Truly, we know a peculiar case when both matrices and are nec- 
essarily positive: that is when A is either in the left or the right gauge. For 
instance, let us assume that A is in the left gauge, then the uniqueness con- 



dition translates into the requirement that the CPT map (2.6) is mixing (see 
If this requirement is satisfied, then automatically: Aq = 1, 
A, but this implies that ie^\e^) = {^^\A) = Tr[A] = 
1 So the denominator in equation (3.20) vanishes. 



appendix A.3 ) 

= Idxd and Cq 
1, telling us that ^ 
simplifying to 



Pe 



J2 ($+|(A., ®A:J...(^,®A:J|A)|si...s,)(ri...r,|. (3.23) 



ri ■ 



■re 



Similarly, the argument can be applied to a right-gauged A, leading to a 



version of (3.23) where correlation boundaries are exchanged. 



3.4.1 Size matters? 

Being able to address infinite states in a formally exact fashion, while em- 
ploying a finite amount of resources, sounds really nice and useful for nu- 
merical issues. Honestly, in many simulation algorithms and architectures, 
the computational times often scale non-trivially with system size L, which 
makes approaches to the thermodynamical limit clumsy trials: expensive and 
imprecise. With Matrix Product States (as well as with other classes of self- 
similar tensor network structures, like CPS, TTN and MERA, we are going 
to describe in the next chapters) the infinite problem is perfectly addressed 
with finite-effort numerics: the amount of calculus scales solely with D, the 
correlation bondlink dimension. 

A good question, now, is whether a Matrix Product Simulation is actually 
capable of describing, with a manageable D, physical states with a good 
precision, even at the thermodynamical limit. We understood that a D- 
bondlink MPS can represent exactly the whole class of finitely correlated 
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states, with entanglement bound by logD; this extends to thermodynamical 
hmit as well. Now consider a ground state of a ID non-critical system: its 
partition entropy ought to satisfy the area-law of entanglement, stating that 
as L grows, Syn scales like L°, going towards a finite value e in the TD-limit. 
Therefore, even the appropriate bondlink dimension D = 2^ stays finite, 
meaning that a thermodynamical MPS description can be definitely made. 
More precisely, in literature there are several classes of quantum states which 
are physically relevant, and an exact (typically optimal) MPS representation 
has been found (see e.g. the review [lOj, where MPS for AKLT, Majumdar- 
Gosh, GHZ, W, and Cluster states are presented). 

On the contrary, this argument suggests us to think that representing 
critical ground states, whose iSvn ~ (clogL)/6 require infinite bondlink D 
for a successful description of the thermodynamical limit. This looks to be 
true, and despite this would be troublesome for numerical implementation, 
the issue of dealing with infinite-bondlink Matrix Product States has been 
studied from analytical perspectives, like in refs. |l3l E] where an equivalence 
with conformal field theory (CFT) has been established. 

To provide more clearance on the relationship between Matrix Product 
representations and criticality in ID, we are going to discuss about correla- 
tions in MPS. 



3.5 Matrix Product States and correlations 

One dimensional quantum systems are quite peculiar: they can manifest 
no quantum phase transition at finite temperature, nor they can exhibit 
long-range order parameters. Yet, the relationship between criticality and 
non-criticality of a system is a matter of utmost importance. As quantum 
entanglement is not a physical observable, the basic way to recognize and 
identify the presence of a quantum phase transition is through scaling laws 
for two-point correlations. 

Let us go back to the thermodynamical limit (TD) MPS, defined ho- 
mogeneously by the matrices Ag, and define the correlation function of two 
(separate) local observables G and G', acting at arbitrary distance i + 1: 

£,+1(0, G') = (0['°] ® 0'[^o+^+i]^ _ (3_24) 

since the system is translationally invariant, io is irrelevant. To workaround 
normalization issues we will just adopt an A in the left gauge, as we did for 
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(3.23), then we end up with 



c:m(0,0') 



(<l>+|Ee E^, Ee'lA) - ($+|Ee|A) ($+|Ee'|A) 
(<l>+|Ee- K-|A)($+|] -Ee^lA) 
(^riK-|A)($+|] K.^f*). 



(3.25) 



Let us now study the expression in (3.25) from an algebraic viewpoint. As we 



previously stated, both left and right (double) correlation boundary vectors, 
!•) and (•!, can be seen as matrices, whose vector version is the so-called 
Liouville representation: 



O = ^aij\i){j\ 



\0) = ^aijij), 



(3.26) 



where, in particular 1 — )■ |$^). Then Ei is a Complete Positivity, Trace 
preserving map when applies to the matricial element on its right, as we 



know from (2.6); similarly, it is a Complete Positive and Unital (i.e. it maps 
the identity into itself) map when applying to the left. 

We begin requiring that the CPT map has mixing property, meaning 
that has a single attraction point, that would be A. CPT maps are always 



contractive, as proven in appendix A. 2 but the relaxation requirement tells 
us also that E°° = |A)($"^| so that we can write 



right X 



(3.27) 



which has a clear physical meaning: the uncorrelated product of expectation 
values is, as should be, equivalent to the operator product at infinite distance. 



We can now exploit the expansion (A.8) for multiple application of a CPT 



mixing map and write the correlator as: 



(r,+i(e,e') = $^A^,pitk^) 



(3.28) 



where Xa>2 are the eigenvalues of Ei other than 1, and for which it holds 
|Aa>2| < 1; while V^^\i) are polynomial functions of degree x with coefficients 
depending on ^, and is the size of the largest Jordan block belonging to 



the generalized eigenspace of a. Two features of (3.28 ) are worthy of remark 



Since every has modulus strictly smaller than 1, €e goes necessarily 
to at £ — )■ oo, regardless of G and G', which is telling us that the state 
manifests no long range order parameter. Notice that this property 
depends strictly on the mixing condition employed for E^. 
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Correlations decay exponentially. In particular it is possible to domi- 
nate €e with decreasing exponentials: precisely, let us order the eigen- 
values {Xa}a so that they decrease in modulus. Then 

i£«(e,e')|_„ (3^23) 



lim 



(|A2|+^) 



0, 



for any e > 0; but if e is chosen small enough (i.e. e < 1 — IA2I) the 



denominator in (3.29) decays exponentially, and thus the numerator 
decays faster. 



It is still possible for two points correlators of MPS, as (|3.28|), to resemble 

a\ very 



power-law decay rates, by playing with several exponentials with |A, 
close to 1, but only for short ranges, and many eigenvalues are required so D 
must be chosen appropriately. When the TD-limit MPS state is investigated 
at long ranges, its ultimate non-critical nature becomes clear, and dominant: 
ruled by the eigenvalues Aq of the Identity transfer matrix E^. 



3.6 Faster expectation values for PBC-MPS 
at large sizes 

When using Matrix Product States as a variational tailored wavefunction 
ansatz for classical simulations of quantum systems, is fundamental that we 
make economy on every computational step of the algorithm. Calculating 
expectation values of observables, and in particular Hamiltonians, is one of 
the numerical ingredients which require most computational effort so it is 
important to optimize its scheme beforehand. 



In section (3.3) we acknowledged that acquiring expectation values in 
periodic Matrix Product States is quite more expensive than in the open 
boundary case, carrying an overall multiplier to the cost (from oc dLD^ 
to oc dLD^). In a paradoxical way, as we approach the thermodynamical limit 
L — )■ 00 the cost drops again to oc dLD^ (plus solving a fixed point equation, 
usually subleading), since boundaries of the TD state are not correlated 
through external channels. Therefore it is natural to think that if the size 
L of the periodic system we are considering is sensibly large, the outcoming 
state shall be close enough to the TD-limit, and thus the system will feel 
little of the periodicity, identified by small amplitudes for finite-size effects. 
A way to exploit this fact to improve the evaluation algorithm was proposed 
in ref . [35] , we now sketch the same idea with a slightly different formulation. 

Assume we want to calculate the composite Transfer Matrix 

P = EW.Egi.....EM, (3.30) 
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where we chose a tensor product observable (^^, ' for simphcity, but the 
following arguments apply to an entangling operator as well. We will also 
state that the whole MPS segment which is calculated through, is 

in the left (or right) gauge, but not necessarily homogeneous. In particular 
if sites 1 ... £ are not the whole system, i < L, it is always possible to take 
the singular part to its complementary, and satisfy the gauge condition. We 
want to achieve but without performing singularly any product E - E which 
costs ~ D^. To this purpose, let us consider the singular value decomposition 
of E^ =, i.e. 

7 

in a formal sense, where WU = = 1, and a|^^ = 6ap era is diagonal and 

positive (cTq, > 0). Now, if the Egj^ were all positive and homogeneous, the 
singular values a would coincide with the eigenvalues and U = V, which leads 
to cTq = A^, where are the eigenvalues of E. But this tells us that the ratio 
between two singular values cxa/cr^ = (Aa/A/j)^ decays as much fast as the 
segment i is long. Telling us that for i long enough very few singular values 
a are relevant before reaching the numerical precision of the calculator. This 
argument holds for positive E, but extends naturally to hermitian matrices, 
and by linearity and continuity it reasonably works for every matrix, for some 
i large enough, and holds even for dishomogeneous matrix products (proven 
numerically in [15]). 

We could equivalently state that the range of E^ has actual dimension p 
smaller than D^, where p is the number of singular values a that are being 
kept, while the other are discarded as they are of equal or smaller order of 
magnitude than computational precision. During simulations, p becomes a 
parameter, and can be kept smaller and smaller as the size i L increases. 
Then the stochastic procedure for calculating E^ goes like this: 

1. Choose a random matrix X, of dimension D^xp. If the random number 
generator is satisfactory, there will be zero probability that one column 
vector will be linearly dependent on the other p — 1 ones, due to the 
fact that those p — 1 may generate a set of zero probability measure 
(absolutely continuous with Lebesgue measure). 

2. Apply and calculate Eq^ ■ . . . Eq^ ■ X = ■ X = Y, obviously starting 
from the right. Every step costs 2dD^p operations. Y is again a x p 
matrix, but the vector columns will span only the range of E^ which we 
will now assume is exactly p-dimensional. As before, chances are that 
all column vectors of Y will be linearly independent, so they will span 
the whole range E^, because dimensions match. 
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3. Orthonormalize the columns of Y, either via a Graham-Schmidt or a 
QR decomposition. We obtain Z = Y ■ T, with T typically triangular, 
containing the whole singular part of Y. The matrix Z is still x p 
dimensioned, and isometric: Z'^Z = 1. Moreover, since the columns of 
Z span Rng(E^), we have Z Z'^ = -PRng(E<) ^'^^ projector over the range. 
But then 

= ^Rng(E^) ■1-' = Z.Z^.I} (3.32) 

4. Apply and calculate Z^ • ■ . . . E^^ = Z^ -f!' = W , from left to right; 
the cost is ~ 2dD^p per step. We are done now, since 

= Z-W (3.33) 

and we have both Z and W matrices. We do not even have to multi- 
ply them together, and instead keep them separated: whenever we will 
have to use E^ as a part of a whole MPS-network contraction, contract- 
ing over the p-dimensioned index space in the middle will be the last 
operation to be performed. 

In conclusion, we can resume these simple steps as follows: 




If we wish to check explicitly the behavior of singular values of E we just 
have to perform a SVD on W, because 



= Z ■ = Z ■ (t/ • aJJ ■ 0) = (Z ■ f/) ■ aJJ ■ 0; (3.35) 
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but {ZU){ZUy = 1, meaning that (3.35) is a singular value decomposition 
for and singular values of a matrix are uniquely defined; so A|^' = A^'^^ 
and we can see how fast they decay for large i. By adopting this process, we 
spent a total amount of elementary computational operations equal to 

#cost ~ 2Lp {2dD^ + SD'^) , (3.36) 

a nice improvement, even because at reasonable lengths {L ~ 50), p can be 
usually chosen orders of magnitude smaller than with practically no loss 



in simulation precision (see figure 3.3) 



3.7 Minimization algorithms 
with periodic MPS 

For open boundary conditions MPS we presented a fast-converging and nu- 
merically manageable algorithm to find the ground state of a generic (short- 



range interacting) Hamiltonian. In section 2.9 we discussed that the essence 
of such algorithm is minimizing one MPS block at a time, keeping the other 
fixed; the basic step is composed by a partial contraction of the MPS net- 
work with the Hamiltonian operator (effective Hamiltonian), followed by a 
c/D^-dimensioned eigenvalue problem {(PD"^ in the two-blocks simultaneous 
minimization case). 

In the periodic boundary case we have more than one naturally available 
path. If the Hamiltonian is translationally invariant, then a good guess would 
be using the set of homogeneous MPS as variational wavef unctions, because 
a translational ground state must exist. This idea would lead to a all-at-once 
minimization of the MPS state, but unfortunately the Lagrangian would not 
be quadratic in the MPS homogeneous block and the problem to solve would 
be way harder than an eigenvalue problem. Moreover, forcing the variational 
state to be translational would let us not identify easily Hamiltonians bearing 
a spontaneous translational symmetry breaking. 

For these reasons, in this section we will instead describe an algorithm 
that not assumes translationality in the variational state (and thus using 
dishomogeneous MPS) and minimizes blocks one at a time to preserve a 
quadratic structure for the Lagrangian [HI |l2]. Then let us start again 
from a nearest neighbor Hamiltonian H = Yle=i -^2' where for comfort we 
regrouped in R2 both one-body and two-body terms, i.e. 

Rf = J29l^^'''^ ®eM + $^/i:ef-^i®e'f. (3.37) 
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Let us assume that we are going to minimize the MPS tensor block as- 



sociated to site i. First, we spht the Hamiltonian as if = -Rf + i?2'^^' + R 



2, 



where R2 contains all the terms of (3.37) that have support in the comple- 
mentary of site i. Then we calculate two composite of transfer matrices, 
which shall be the ingredients of our Lagrangian functional, namely 



The transfer matrix Ei from site £ + 2 to 
i.e. El = Ef^l ■ . . . ■ Ef-^l 



2 of the Identity operator. 



The transfer matrix E2 of R2 from site £ + 1 to £ — 1. Even though 
R2 is not a separable operator it is possible, with some engineering, to 



calculate E2 spending quite the same computational cost (3.36), apart 
a non-scaling prefactor. 



Thanks to the technique (3.34), acquiring these transfer matrices is efficient; 



actually we prefer to store in memory Zi, Wi and Z2, W2 (where E^ = 
ZaWa), since p ^ so keeping 4{px D"^) matrix elements is less expensive 
than 2 [D"^ x D"^). 

Then the Lagrangian for reads: C{A^^\A*^^^) = {^\H\-^) - e{^\^) = 
((AM|7{-£A/'|AM)) = 



c = 




£ 




(3.38) 



M = Ao (8) tdxd is the effective square-norm operator, where A/q is given by 





(3.39) 
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while the effective Hamiltonian "H is obtained as follows 




(3.40) 

We immediately see a difference from the corresponding Lagrangian in the 



OBC case (2.86): the effective (square) norm operator is no longer the iden- 
tity Af ^ 1. If we recall correctly in the OBC case it was a property strictly 
depending on the choice of a gauge condition for the other MPS blocks. But 
when the MPS design is Periodic, in general there is no trick with gauge 
transformations in order to map A/" into 1. This also means that to find 
the optimal A^^'^ one has to solve a generalized eigenvalue problem, instead 



of a simple one as in (2.88). Precisely, the Euler-Lagrange equation of our 
problem is: 

a£(AM,A*M; 



|o)) 



(3.41) 



whose solution with minimal e is the optimal one, since e is actually the 
energy: e = ((^M |H|v4M))/((AM lA^I^M)) = After we found 

the solution, we put the optimized A^^l into the MPS and choose another site 
for the variational paradigm. Then repeat and sweep until convergence is 
achieved. 



3.7.1 Stabilizing the generalized eigenproblem 



Dealing with the generalized eigenproblem (3.41) is no small trouble. Even 
with the most advanced linear algebra techniques, numerical costs are much 
greater than those required for addressing standard eigenproblems (for com- 
patible dimensions). For the latter, eigenvector related to minima or maxima 
of the spectrum are often found quickly thanks power method-inspired proce- 
dures. Several algorithms for the generalized eigenproblem, like the Jacobi- 
Davidson, are also based on power-method principles, but they should be 
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# sweeps 



10 



Figure 3.2: Offset E* from the ground state energy of a variational periodic 
MPS a function of the number of algorithm sweeps. The model considered 
is a ID spin-| XXY ring, with anisotropy A = 0.5 and twisted periodic 
boundary conditions (0 is the twisting phase). Here a bondlink D = 18 was 
used. The graphic was kindly contributed by D. Rossini [I2] . 



wielded with care. These protocols are very efficient when the inverse matrix 
oi J\f is available; and, when it is not, converge faster the more A/" is easy to 
invert. Obviously, an operator A/" is not suitable for inversion when there are 
eigenvalues which are much smaller (closer to zero) than other ones, because 
numerical methods perceive the relative eigenspaces as if they were a kernel. 
In this sense, we can relate the 'fast-invertibility' of a (positive) matrix with 
the requirement that the relative spreading of its eigenvalues is small, i.e. 
(Amax — Amin)/Amin "C 1, and also mcaus that Af is somehow 'close' to the 
identity, as the multiples of 1 are the only positive operators having relative 
spread zero. 

Now, can we perform some gauge transformation that takes A/" as close 
as possible to the identity? We stated that there exist no general solution 
to this question, although, for large system sizes £ we can argue that the 
system shall manifest small finite-size effects, and an affinity with the OBC 
version should be met. In section [3l6] we saw that the largest eigenvalues of 



a composite transfer matrix decay more fast the larger is L. In particular, if 
every Ei was mixing-CPT, we could write 



E{ \Al){Ar\ + e' 



(3.42) 



where \\0i\\ is bounded regardless from £, and e is somehow related to the 
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largest (in modulus) eigenvalue of E smaller than 1: e ~ |Aq,| < 1. It is 



clear that for i oo the second term in (3.42) vanishes, but also for i finite 
but large Oi is just a small perturbation to the first, leading, term. A^, 
and Aij, read as matrices, are positive thanks to CPT condition. 
^ Relying on this concept we^perform the following operations, consider 
Eg = eJ+^I ■ . . . ■ eJ"^' = eY+^^ ■ El ■ Ejf"^^; if L is large and all the fixed MPS- 



blocks are in the left gauge, then E3 is written as (3.42) and in particular 
(A^l = ($"*"!, Then, since A^ is positive, we can write (via SVD, for example) 
A^, = XX"!", or equivalently |Ai) = (X ® X*)|$"'"). Now we can perform a 
gauge transformation upon A^^"^! so that it adsorbs the X operator: Ai^^^' — )■ 
A^s ■ X. After the transformation, we are left with 

Eg ~ !$+)($+! +£^0;, (3.43) 

which, in turn, makes the effective square-norm operator to read as follows: 

A/"- l + e^TVj^^, (3.44) 

whose relative spread of the eigenvalues scales with e^, where e < 1. We 
ended up with an effective normalization operator which is actually the iden- 
tity apart a small perturbation, which decays exponentially with the length 
L. When we apply the Jacob i-Davidson method in this framework, we find 
the generalized eigenproblem solution much faster than if we do it naively, 
as proven in numerical simulations [HI |12] ■ 



3.8 MPS and Tensor Networks 

With this last consideration, we conclude our discussion on how Matrix Prod- 
uct States (either in their open or periodic boundary contexts) relate to 
simulation paradigms as variational tailored wavefunctions, with surprising 
efficiency even at high precision calculus, and wide manipulation features 
that let them overpower the old-fashioned DMRG design. 

As you could imagine, though, their application does not limit to numer- 
ical settings. MPS are kept in great regard even for analytical calculation 
purposes. The possibility of building interesting parent Hamiltonians for any 
Matrix Product State [lO], and their continuous-space enhancement to de- 
scribe a finitely-correlated quantum field theory [12], HH], are just two of the 
many developments achieved in the last ten years. 

More than anything, MPS have been the main inspiration that led physi- 
cists to investigate thoroughly in the simulation capabilities of Tensor Net- 
works in general. PEPS, TTN, MERA, and other well-known Tensor Net- 
works designs, actually gathered interest only after the MPS representation 
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Figure 3.3: Convergence rates of the energy offsets with sweeps, using the 



same ID periodic XXY model of figure 3.2 , for various transfer matrices trun- 



cation parameter p, as explained in section 3.6 and equation (3.34). Notice 
that p = 30 achieves an excellent precision, even though the original double- 
bondlink had dimension = 324, thus a factor 10 in algorithm speed. The 
graphic was kindly contributed by D. Rossini |12]. 



of DMRG states was fully understood. In the following chapters we will try to 
understand how such Tensor Network-based structures work; how the ideas 
behind MPS technique can be redesigned in order, for instance, to extend the 
success of this method even in settings where standard MPS/DMRG fails, 
like higher-dimensionality systems. 
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Chapter 4 



General features of Tensor 
Networks 

Every many-body /multipartite quantum state, when described as amplitude 
components over a separable product basis, is uniquely represented by a 
complex Tensor. However, the number of identifiers, of descriptors, one has 
to express to locate specifically that state, within the whole manifold of 
system states, is huge: in principle it scales exponentially with the number 
of elementary constituents of the system, regardless of their nature. Tensor 
networks (TN) are the trial to express the same state with a number of 
numerical descriptors (be they variational or parametric) which is small, 
meaning that they must scale nicely with the system size L, and that lead to 
the original wide amount of descriptors by means of simple linear algebraic 
operations. It is really not fundamental whether the desired analytical state is 
reproduced exactly as much is important to recover the real-physics features 
the state should exhibit. 

It is really impressive to acknowledge how much interesting physics can 
be generated for an approach that sounds so naive, even in contexts where 
other analytical or numerical methods are totally clueless. In this chapter, 
inspired by what wc learned about MPS (the father-archetype of TN) we will 
try to understand what common properties and features the Tensor Network 
architectures share, and also present some useful examples and comparisons. 

4.1 Definition of Tensor Network state 

Let us start from a generic multipartite system "H®^ = (3)t=i ^> where "H 
is the Hilbert representation of the elementary degree of freedom; here we 
are taking the constituents to bear a homogeneous representation, as it is 
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the usual physical setting, but this is not really a requirement of the TN 
description. An orthonormal basis for the single degree of freedom must be 
chosen as canonical one, which we shall represent as |s); then it is standard 
procedure to expand any given state of T-L over the product canonical 
basis: 

d 

|^)= Yl '^^--^ \siS2...sl) (4.1) 

S1...S£=1 

where d is the dimension of the elementary Ti, L the total number of con- 
stituents, and Ts^...sL is a complex tensor (i.e. a multi-indexed collection 
of complex numbers), with L indices, each one allowing d different values. 
The normalization condition reads ^j^^}^ '^i...sl'7^*...s^ = 1- As long as we 



ordered completely the L degrees of freedom, (4.1) is meaningful for every 
dimensionality of the physical system, and for every nature (spin / bosonic / 
fermionic) of the constituents as well. In fact, the complete ordering allows 



us to write always a second-quantization version of (4.1), as follows 



l^)= E r.,...., (cI)-(4)-...(cl)^^|(]). (4.2) 

Sl...SL = l 

where c]^ are either Bose or Fermi operators, and in each case they satisfy the 
proper commutation or anticommutation relations. Spin-orbitals a are now 
the elementary components of the system, their on-site filling being the local 



canonical basis. Expression (4.2) represents the more general many-body 
state, and Ts^...sl contains all its physics and information. 

At the same time T is a huge array, with d^ elements, hard to manipulate 
in every sense. But assume that the state is such that T can be obtained, 
via contracting over an additional index q, from a pair of tensors T' and T", 
like 

T=^T' T" . (4.3) 

' / ^ ' si...si>,q ' Si+i...si^,q- \^-'-'J 

9 

Then we would need a number of descriptors equal to D{d^ + d^~^) ~ Dd^^'^, 
with a lot less information needed when the number D of allowed values for 
index q (or index q dimension) is smaller than d^/^. The same argument can 
be applied again and refrained for T' or T"; every time we split a tensor into 
smaller tensors partially contracted together. Every time, there is a chance 
(and typically happens) that we lose description capacity, meaning that the 
resulting set of states allowing the new decomposition is often than before. 
But this is not an issue as long as the states we are cutting out of our ansatz 
are those which are not physically relevant, and we keep those that contain 
the true physics of the problem we want to study. 
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This is the central point of the Tensor Network paradigm, inspired to 
MPS. We write T as product of muhiple tensorial objects, where L indices 
are left open, they are the d-dimensioned physical indices, while an arbitrary 
number of 'fictitious-space' indices {q} (of arbitrary dimension Dq) are con- 
tracted. Like for MPS, for a given scheme of contraction, which from now on 
we will call Network, there is always a choice of virtual links Dg large enough 
so the whole is described. But that would be a waste of effort, since 
it is very unlikely that a physical state would require that very amount of 
information shared among tensors in the network. Which yields a definition 
of TN-state that takes into account parametric bounds to our description in 
terms of (i) the number of tensors, (ii) the number of indices per tensor, and 
{Hi) allowed values per non-physical index: 

Definition of Tensor Network - A multipartite state is a Tensor 
Network state {V,k,D}, with maximal tensor number V, correlation link- 
number k, and link dimension D if 

• It exists a decomposition of Tsi...sl as a contracted product of tensors: 

q- _ \ ^ \ M I I T^M'J«'(l.")--'?«'(''n,n) 1 

lsi...SL — / ^ ■ ■ ■ / ^ I II | ) l^-^J 



where V is the total number of elementary tensors = nodes in the 
network, B is the total number of contracted indices = connected links 
in the network. Every appears once as tensor index in the expanded 




expression (4.4), and every appears twice 



• The number of tensors in the decomposition is bounded by V < V, 
the bondlink dimensions are bounded by < D for every a G {1..B}, 
and the total amount of indices of a tensor is bound by a„ + 6„ < A; for 
every n G 

As an immediate consequence, the overall number of complex value de- 
scriptors (i.e. variational parameters) required for such representation is 

^elements < V ( sup{L), d})'' . (4.5) 
Although this is a well-formulated definition, there is no doubt that the for- 



malism of (4.4) is cumbersome and confusing. For most purposes involving 
Tensor Network states is actually preferable to involve a diagrammatic rep- 
resentation, much similar to the one we adopted for Matrix Product States. 
In these diagrams. Tensor Network states are represented as graphs: Tensors 
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being vertices, indices being links, either left open if they are physical indices, 
or connected if they are fictitious. Let us consider the following example: 





L 

(4.6) 

This network is made out of = 8 tensors, with maximal correlation number 
k = 6; the analytical tailored expression of the corresponding (spin) state is 
given by 

d Dfi<D 

|v]/\ = y[2]'7i<?2g3 2^[3]g2g4g5 2^[4]i33g4g6 

{sa}=i {g/3}=i 



which is messy, and not immediate as (4.6) although they represent the same 
parametric set of states. 

Most classes of Tensor Network commonly considered in literature, are 
scalable with system size. We intend that the network is obtained by re- 
peating some fixed local pattern of vertices contraction to build a self-similar 
structure, so it can be adjusted to fit any L bu just adding new tensors ac- 
cording to the same pattern as before. When doing so, it is important that 
k and D can be kept fixed, and the number V of tensors (recall that V is 
proportional to the number of variational parameters) scales nicely, e.g. lin- 
early, with L. This was the case of Matrix Product States, where precisely 
V = L. D is also occasionally referred to as refinement parameter [47J, as 
its value directly influences the capability of the TN ansatz. 



4.2 Entanglement of Tensor Network states 

Quantum entanglement is the primary responsible for argumenting that Ten- 
sor Network are actually a good technique to describe physical multipartite 
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states. Ground states are characterized by small correlations, compared to 
a generic random state in the Hilbert. So probing them with trial wave- 
functions that admit a simple description, and yet capable to reproduce just 
the needed amount of entanglement seems a suitable choice. Tensor Net- 
work states have the outstanding feature that their entanglement is perfectly 
controlled by the network topology itself, as we will show in this section. 

Entanglement bounds of a Tensor Network state - Assume 



is a state which allows a Tensor Network representation, as in (4.4). Let 
us choose any partition of the physical sites sj into two disjoint subsets s^^^ 

r D] 

and Sj . Then the Von Neumann entanglement entropy associated to this 
partition satisfies the following inequality: 

5vN(pt^l) < min I VlogDj (4.8) 

where the minimum is taken over all the partitions of the network graph into 
two subgraphs, with the condition that the first subgraph embeds all the A 
sites and the other one the B sites. The are the bondlink dimensions of 
the links we need to break in order to disconnect the two subgraphs. 

To make an example, let us consider again a Tensor Network state like 



(4.6) (4.7), and assume that we are to estimate the entanglement shared 
between the six leftmost sites A = {1..6} and the five rightmost sites B = 
{7..11}. Then identify all the possible ways to part the network into two 
subnetworks, respectively containing A and B. Three smart choices are given 
by: 





(4.9) 

Finally detect all the network bondlinks a we should break to separate the 
two subgraphs (violet and green), and sum their logD^ to obtain a bound 



on the entanglement. In particular, the three graph partitions picted in (4.9) 
tell us that 



Slf^""^ < imn{\og{DsD,D,),\og{D2D,),hg{DjDs)} (4.10) 



and, since any Da < D, we conclude that 5vn < D^- The reader can easily 
check that no other graph partition would lead to a tighter bound. 
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Let us prove the statement (4.8), by adopting an argument very similar 
to the valence bond picture for MPS. Consider any network graph partition 
into two subgraph, in accordance with the lattice sites partition A ^ B, a.s 
above. Then, let us write a starting virtual state 

1*^) = (g) \KJ where |$+) = ^ E l^^)" (^.ll) 

Each maximally state |$J ) contributes to the entanglement of sepa- 
rately, since they lie in different degrees of freedom, and each one contributing 
with log Da- But now we can find a linear mapping taking the state |$^) 
into the original l^') 

|M/) = (Affl®Mf)|$+) (4.12) 
the mapping m|^' is given by the contraction of the subgraph on A, and 



therefore is linear, and so is M[ telling us that ( |4.12D can be viewed as 
a quantum transformation, not necessarily invertible. But as it is a tensor 
product, it is also local, and thus can only degrade entanglement, not enhance 
it. Therefore, the entanglement of |\E') must be less than that of |$^) which 
is exactly log Da- The same argument can be repeated for any subgraph 
partition, thus concluding the proof. 



4.3 Operators and link exchange-statistics 



Since we now have a diagrammatic representation for linearly-connected tai- 
lored variational wavefunctions, we want now to exploit this idea also to 
include the action of operators. Especially operators that act locally or at 
short-ranges appear as new tensorial pieces to add to the network structure, 
connecting to those physical links that were left open in the TN-state design: 
for instance, a three-site operator (acting on sites number 2,3 and 4) applied 



to (4.9) reads 



e 



2,3,4 



I*) = 




(4.13) 
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As was pointed out in refs. [HI HHl |50], attaining this representation is un- 
doubtedly trivial for a spin or boson system, but we should handle the issue 
with care when our Tensor Network state is describing a system of fermions. 



Indeed, when we introduced the second quantization version (4.2) of Tensor 
decomposition, ordering the sites was a crucial point. It is clear that chang- 
ing the ordering of those sites would not only rearrange the components of 
T but also change some signs appropriately due to Fermi statistic. We are 
about to show that this feature can be embedded as an inherent property 
of the network links, that manifests under crossing (exchanging) of the links 
themselves. 

For instance, we will start from a two-site operator Qix, we will also 
assume for simplicity that this operator preserves parity, i.e. 

9i_2 = ao 1 + «! C1C2 + a2 c\c2 + 03 cicl+ 

+ a4 c\c\ + as c\ci + clc2 + c\ciclc2 (4.14) 

requiring parity preservation has the advantage that (quasi-) local opera- 
tors preserve support locality when representing fermions as spins. Now, 
assume the tensorial representation of Bi,2 is known and available, encoded 
through descriptors {ao . . . 04}. We wonder: how do we express the action of 
9i^3 instead? The standard procedure is swapping the second and third (or 
equivalently, first and second) sites both before and after performing 0i^2; as 



S 



( 



e 





e 



) 



(4.15) 

According to anticommutation rules, it is during the exchanging process that 
the fermionic statistic should emerge. Indeed the operator S in equation 



(4.15) is the standard quantum information Swap gate only for spin and 
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bosonic Tensor Network. For fermionic Tensor Networks it reads instead: 




/I \ 

1 

1 

V -1/ 

(4.16) 

where the antisymmetrization sign —1 appears in correspondence to cjcglfi) = 
— C2cJ|f2), the doubly occupied canonical state, all the other 2-sites canonical 
vectors, namely c\\Q) and C2|fi), being either zero or one particle states, 
feel no difference from the spin/Bose setting. 

In conclusion we acknowledged that Tensor Networks ansatze can be suc- 
cessfully applied to fermionic systems with no basic difference in efficiency or 
computational costs. The only issue we need to take care of, is the exchange 
of network links within diagrams, that can be rearranged from the original 



ordering by means of (4.16). We could state that the network links satisfy 



an exchange statistics themselves, that identifies the nature of particles the 
TN-picture is representing. 

As an additional remark, let us point out that it is straightforward to 
generalize all these arguments to include abelian anyon (like those of frac- 
tional quantum Hall effect) algebras as well. Indeed, let us assume that we 
are representing, with our trial Tensor Network state, a 2D system, which 
is where anyons arise. Then if anyonic operators undergo the exchange rule 
then the the bondlink exchange statistics of the TN repre- 



C1C2 



sentation (4.16) is replaced by 



S: 



anyonic 



/ 1 



V 





1 



(4.17) 



The problem of extending such formulation to include also non-abelian anyons 
was discussed in ref. |51]- A similar scheme to deal with fermionic statistics 
is by using directed-link Tensor Network designs introduced recently, called 
Fermi operator circuits (FOC) [50]. 



4.4 Gauge group of Tensor Network states 



Inspired by our discussion of section 2^ where we defined a group of MPS 
transformations under which the physical state is invariant, we want to ex- 
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tend this paradigm to any given Tensor Network geometry. As for Matrix 
Product States, a gauge transformations group provide a clever way to ma- 
nipulate analytically or numerically tensors within the network structure. 
The purposes of exploiting this group are many, for example, to perform 
faster contractions. 



The generalization of (2.44) to an arbitrary TN-state is natural. Consider, 
for instance, a shared (connected at both sides) bondlink z/, corresponding to 
a virtual degree of freedom, of dimension D,^. Tensors T^"' and T^^l are the 
two nodes in the network sharing link u. Then we can define the composite 
tensor emerging from the contraction of the two 

j,[a-n-/3]{qa,gp} _ Sr-^ rp[a]{qa},q^ rp[P]{qp} ,qi' j^gN 

We are also assuming that no homogeneity or geometric constraint is re- 
quested for the TN to hold, so that the state |\E') actually depends on T'"! 
and T^^'^ only through T'"^'^', i.e. even if each of the former two changes but 
the latter is unaltered, then l^') is also unchanged. 

Now let us choose an isomorphism X on C^", i.e. an invertible linear 
application on the D,^-dimensional complex vector space: XX~^ = X~^X = 
^D^xD^- Then it is clear that 

rp[a^l3]{qa,qi3} _ rp[a]{qa},ri y rpll3]{qi3},r3 

— - (4,19) 



EyyH{qa},qu ^rW]{Q|3},q^ 



91^ = 1 



where W^"^ is the contraction of T'"' with X, while W^^^ is obtained by linking 
X-i to as in the following diagram 




(4.20) 
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Since T^°'^^^ is invariant under this transformation, is invariant as well. 
Rephrasing (4.19) in other words, we could say that, as u is a fictitious, 
external degree of freedom over which we are performing a (partial) trace, 
the physical system (set of real degrees of freedom) is insensitive to any lo- 
cal invertible transformation X acting upon z/, which is what (4.19) does. 
There is no need to say that transformation (4.20) has a natural group struc- 
ture, arising from the fact that the set of Isomorphisms {X} is closed under 
composition. 

This very argument can be applied and repeated for every closed (doubly 
connected, i.e. non-physical) link in the network structure. Indeed, we asso- 
ciate an invertible matrix X,^ and a direction to every virtual TN-bondlink 
and transform tensors T^"! — )■ W^°''^ according to 



'a]{q},{p} 
{s} 



n 



X, 



T, 



'a]{v},{w} 



(4.21) 



{v},{w} 



where k^^ (resp. k"^^) is the number of links whose chosen direction is in- 
coming to (outcoming from) node a. Precisely, we are contracting T'"! to X^ 
at link z/, if the direction given to u is outcoming from node a, and to X~^ 
otherwise. As stated, no transformation is allowed on the physical links, i.e.: 




It is worth mentioning that the group we found, which coincides with the 
Cartesian sum Q = ^"^0(0^") (with n the total number of closed links), 
is the most general state-invariant transformation, provided that the TN- 
geometry and the bondlink dimensions Dy are fixed. 
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4.5 No closed loop ^ efficient contraction 

As we learned how to match a quantum operator with the Tensor Network 
representation of a state, we need to work out how to achieve expectation val- 
ues, to get physical information on the state itself. For MPS we argumented 
the importance of a contraction algorithm which is numerically efficient: as 
we typically adopt energy as a benchmark for simulation convergence to- 
wards the ground state, it is likely that we will need to evaluate 
many times. 

An optimal contraction scheme must be modeled on the topology of Ten- 
sor Network we are considering, unfortunately, even when such scheme is 
optimized, its efficiency could still be unsatisfactory for practical purposes, 
and that obviously depends on the network design. Precisely, even if we have 
a L-scalable Tensor Network structure {V{L),k,D} with a number of vari- 
ational parameters V proportional to L, the optimal contraction cost could 
still scale exponentially with L. This is, in particular, the unfortunate case 



of Product Entangled Pair states, which we will introduce in section 4.7 

However, we can identify a sub-class of Tensor Networks which is still 
quite general and whose contraction efficiency is guaranteed. We are talking 
about Tensor Networks that lack closed loops in the graph structure. Ac- 
quiring an expectation value out of a Tensor Network state {V, k, D} without 
close loops has a computational cost bound by 

#cost~2VD2^ (4.23) 

which scales linearly with size ii V cc L. Here we present an example of the 
comparison between TN-structures without or with closed loops: 



where orange dashed lines highlight the minimal loops. 



(4.24) 
The contraction 



scheme that leads to (4.23) is definitely intuitive and often close to be the 
optimal one for a given network. Let us sketch it briefly. 
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Let 00 = <S)^e be the separable product operator whose expectation 
values (\l/|0,g,|\E') we are interested in. Then, for every node fi in the network, 
we define a new 'double'-tensor 



E 



Yl{rvU)\QvU)\svU))- (4.25) 



where n is the number of physical links connected to node /i, which corre- 
spond to sites f (1) . . . v{n), while m is the number of virtual links; of course 
m + n < k. The tensor Wt^^ is obtained by pairing together tensors T''^' and 
its complex conjugate T* ^\ while contracting every physical link though the 
action of B^; notice that W^^^ has correlation number 2m < 2k. Obtaining 
every W^^^ has a computational cost of at most D'^''. 

Now, it is easy to see that the >V''^1 form again a Tensor Network, which is 
identical to the original one except for having no open (physical) links, which 
vanished. So it contracts completely into a real number, equal to (\E'|60|\1'). 
The resulting bondlink dimension is the original one squared: D' = D^. We 
can now exploit that this new network has no closed loop, then we start 
contracting from terminal nodes, i.e. nodes that have correlation number 1 
(a finite graph without close loops has always two terminal nodes at least). 
Contracting a terminal node has a cost in elementary operations that scales 
like D^'^, and i have to repeat this operation ~ V times. Thus the overall 
cost of the full contraction is ~ VD"^^, and by adding the expense to obtain 
the double-tensors >V''^1 we recover (4.23). 

It is important to notice that such procedure was explained disregarding 
the ordering of sites, so it works for Tensor Network states of spins and 
bosons, but it is not so trivial for fermions, since link-exchange operators are 
not separable operators. Nevertheless, it was pointed out in [52j that this very 
scheme can be generalized to the Fermi case, with equally-scaling efficiency 
rates, by adopting Z2-symmetric tensors in the network, with practically no 
loss of generality (see appendix [B]) . 



4.5.1 Peripheral gauge 

We gain a remarkable computational speed-up by manipulating with gauges 
a Tensor Network state without closed loops in its structure. Assume, for 
instance, that we want to acquire the expectation value of an operator 6 
having support on a restricted set of sites {at}. Now, we define the nucleus 11 
as the smallest connected subgraph of the network containing all the external 
links {at}. Then it always exist a gauge for which the number of operations 
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to obtain (6) is 

#cost ~ 2#(n)D2fc, (4.26) 

instead of 21^0"^^, where #(n) is the number of tensors in subgraph 11. This 
is definitely an interesting improvement, especially when the support of is 
very local with respect to the network geometry, i.e. when #(11) ^ V. 



This is the same principle that, in OBC-MPS context, lead us to (2.71); 
and the definition is totally similar. Let us briefly sketch the idea on how 
to convert the tensors not belonging to 11 into the peripheral gauge, that 
eliminates them automatically (with no need of numerical calculus) when 
acquiring (B). First, let us associate to every tensor in II'^, the complemen- 
tary subgraph of 11, a distance: the graph distance to the nucleus 11. Then, 
starting to the tensors having highest distance, we perform Singular Value 
Decompositions 

7 

where q (resp. p) indices are related to links that decrease (increase) distance, 
i.e. that lead towards (far from) the nucleus 11, and s are physical link 
indices. Then we damp all the singular part U ■ X of every tensor in 11^ into 
the only linked tensor with shorter distance (this is a gauge transformation 
mapping T ^ V). This is done recursively, until all the singular part has 
been embedded into 11. And since all the periphery 11'^ is isometric, it cancels 
out analytically when contractions with (^'l are made. 

It is easy to see that only for non-closed loop Tensor Network a peripheral 
gauge can be defined regardless of the support of 9. if 11'^ contained a loop 
then there is no gauge that contracts automatically the branch. This is what 
happened exactly for periodic MPS, whose single loop is sufficient to break 
down the definition of a peripheral gauge. 



4.6 Energy minimization techniques 

Performing the variational step towards the ground state is often the bot- 
tleneck of any simulation algorithm. In Tensor Network ansatze, a good 
procedure of minimization makes the practical difference between better and 
worse network designs. However, there are some common ideas involving 
variational algorithms around TN: first, the globular structure of a state al- 
lows us to treat various tensors in the network as independent variables, and 
therefore we are highly encouraged to minimize one (or few) of them at a 
time, while keeping the other ones fixed. This, of course, is possible as long 
as we do not require some homogeneity constraint among tensors: we might. 
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for instance, wish to fix some tensors to be identical, e.g. due to the presence 
of a global symmetry (like homogeneous periodic MPS under Translational 
invariance); in this case the minimization of various network nodes would be 
simultaneous. Then we should treat the two cases separately. 

No homogeneity constraint - if we are free to adjust every tensor 
independently of the others, then addressing the problem is simple because 
the Lagrangian functional is always quadratic on the single tensor we choose 
to variate C = (\E'|iJ|\E') — £:(\E'|\E'). In practice, we could contract the whole 
network except for the tensor T of interest, which appears twice in the ex- 
pression, once as T and once as T*, similarly to the PBC-MPS case. Then 
£(T,T*) = {{T\H-eU\T)) = 



c = 




e 




(4.28) 

where the effective Hamiltonian "H is hermitian, and the effective square norm 
is a positive operator A^. This is due to the fact that partially contracting 
the Tensor Network around an operator O is formally equivalent to 



tM 



(4.29) 



but this is a completely positive map (see appendix |Aj), which preserves 
operator positivity and hermiticity. Then we are allowed to address the La- 
grangian problem (4.28) as a generalized eigenvalue problem "HIT) = eA/'|r), 
whose minimal e solution defines the optimal T. 

With homogeneity constraint - when two or more tensors in the 
network are chosen a priori to be identical, but not fixed, then the problem 
becomes more complicated as the Lagrangian reads 



£(T, T 




((Tl, [U - sN] 




\T)), 



(4.30) 
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where w is the number of times the same variational tensor is repeated in 



the network structure. The problem (4.30) can be approached with various 



numerical methods those that gathered most interest in literature are 
the following 

• Gradient methods - this idea involves guessing random variations to the 
original tensor by using the Lagrangian gradient dC{T, T*)/dT (or the 
conjugate gradient, which converges usually faster) as favored direction 
of random walk. 

• Linearized problem - here we try to treat £ as a linear functional 
C = {{Q(^T,T*)\T)) whose solution is found immediately by polar de- 
composition. Doing so is very cheap and fast-paced but convergence is 
not guaranteed in general. 

• Imaginary time evolution - Another way of proceeding is starting from 
a random state and cooling down the system by applying e~^^ to it, 
for /3 — )■ cxo the state thermalizes at zero temperature, i.e. is the ground 
state. But writing an unperturbative exponential of the Hamiltonian 
operator H is usually hard task, and not always suitable to match the 
Tensor Network geometry. 

This mostly concludes the general features of Tensor Network states that 
we wanted to discuss in this chapter. Before in-depth investigating the pe- 
culiar properties of the famous class of Hierarchical Tensor Networks (Tree 
networks, Multiscale entanglement renormalization ansatz) we focus on some 
other archetypal examples that are worth mentioning. 



4.7 Example I: 

Product Entangled Pair States (PEPS) 

Product Entangled Pair states are the natural generalization at dimension- 
ality > 1 of Matrix Product States. The idea behind their formulation is to 
adapt the valence bond picture of MPS to 2D or 3D systems. We discussed 



in section 2.2 that an MPS can be seen as an on-site transformation ap- 
plied to a starting, virtual state made by several maximally entangled pairs 
l^"*") = (with D is arbitrary): every pair belonging to a physical 

bond in the lattice. 

It is clear that this very construction is meaningful for any dimensionality, 
and for every lattice we can think of as long as nearest neighboring sites are 
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well-defined. For instance, when we apply this framework on a square lattice 
we obtain: 




(4.31) 

(picture taken from MPQ-Garching) where red dots are physical sites. Every 
site has four bonds, and every bond carries an entangled pair. Then a linear 
transformation maps the D^-dimensioned fictitious degree of freedom into 
the d-dimensioned physical one, that defines uniquely the PEPS. 

Since PEPS are fully parametrized by such local transformation, they 
are actually a Tensor Network {L,k,D}, with a number of tensors equal 
to the total number of sites L, correlation number equal to the number of 
neighbors per site plus 1 for the physical link (e.g. k = 5 for a square lattice), 
and arbitrary bondlink dimension D. Their TN representation reads 




(4.32) 



for a square lattice, where vertical links are physical indices. Red boxes are 
tensors ^if, 92,93, 94) where jj is the site, s the local canonical level index, and 
gi . . . g4 the four correlation indices. 

The reason why PEPS gathered so much interest in the latest years fi5[ 
153] , is that they obey the exact entanglement area law for their geometric 
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dimensionality; e.g. a 2D PEPS satisfy a 2D area law of entanglement. It is 
quite easy to check it, by using the TN-entanglement bound rules we derived 



in section 4.2 Consider, for instance, a connected subset of sites in the square 
lattice, and the entanglement entropy related to the reduced density matrix: 
we mentioned that this entanglement is bound by the minimal number of 
links one has to cut to separate the graph into two subgraphs (times logZ^). 
Then one sees immediately that this minimal number coincides with the 
perimeter of the subset of sites. But that is the 2D area law: for a region 
scaling like a surface ~ i"^, the entanglement, 

SY^^ilogD, (4.33) 

scales with the characteristic lenghtscale i of the region. This suggests that 
PEPS are optimal tools to approximate ground states of short-ranged Hamil- 
tonians, even because, similarly to their one-dimensional cousins MPS, they 
also yield a completeness theorem. 

PEPS have also a huge drawback when compared to MPS, which is their 
computational complexity [54j. Namely, if we want to contract exactly a 
PEPS, e.g. for acquiring an expectation value, the number of elementary 
operators needed scales exponentially with the system size, no matter the 
contraction scheme. This is due to the presence of an extremely high amount 
of loops in the network, not an easy trouble to work around. Of course one 
can deal with the problem approximately, by, for example, renormalizing the 
width of the resulting composite bondlink during contraction; still, we must 
be careful that the approximation we are adopting should not break down 
the area law principle. To this goal, several proposals have been introduced 
in literature [161 EH] • 



4.8 Example II: 

Correlator Product States (CPS) 

Correlator product states [23, ESI EZl EH] are a simple, but clever, way to 
adjust correlations in many body states, by means of a whole abelian al- 
gebra made of ranged canonical weight-factor operators. They have been 
acknowledged and used long since; in some sense, they can be regarded as 
the completely variational generalization of Jastrow factors. 

They are called 'states', but it would be more suitable to regard correla- 
tor product elements as operators. Indeed they hardly stand alone as a mere 
Tensor Network structure, because as they must be finite-support and com- 
muting, they are forbidden to establish global symmetries, even the simplest 
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ones like particle number conservation. On the positive side, they excel in 
preserving symmetries already present in the state prior to their application 
as operators, so they are often used in conjunction to a starting ansatz: 



1^ 



CPS/ 



$0 



|si . . . Sl). 



(4.34) 



Sl. 



In this representation, |$°) = XI ki • • • -sl) is exactly the starting 
ansatz state, capable of controlling desired symmetries; e.g. nonzero ^ 
components should be only those with the right particle number. It is clear 
that correlator product factors Ct"' can not break such symmetry, since sec- 
tors with zero amplitude will stay zero amplitude. Nevertheless, C^"' can 
build correlations, which is especially useful when l^'') is a somehow un- 
correlated trial state, like a Slater determinant (in particular, see the MPS 



representation for Slater determinants in section 2.11). 



Notice two facts involving (4.34): first, the action of correlator product 



factors C'"] can be applied in any ordering, because they are all diagonal in 
the canonical basis, so they must commute as they share a basis of eigenvec- 



tors altogether. Secondly, (4.34) is clearly a Tensor Network, since all the 
relations between the various C'"' and $° are linear; so we can use arguments 
treated in this chapter to study CPS (e.g. to guess entanglement bounds). 
You could argue that such picture does not match exactly the definition of 



TN we gave in (4.4), since the same index is shared by more tensors than 



merely two, but this is easily worked around by adding triple Kronecker delta 
nodes in the network Slj^i^ = 6ij6j^k- 

Precisely, consider a two-site diagonal operator 6, defined as 

(4.35) 

then we can exploit its diagonal nature to decompose its operator represen- 
tation even in the diagrammatic scheme 



ri 



r2 



e 



Sl 



S2 



(Jll^2] 



(4.36) 
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where the black nodes are 5^^^ tensors; so the network structure of (4.34) is 
well-defined. Of course, correlation capabilities of CPS depend strongly on 
the specifics of chosen C'"' factors: 

• Involved sites (k) - any number of sites can be chosen to appear as 
indices in a single C'"! factor. The amount of variational descriptors 
and thus the speed of simulation algorithm is extremely sensitive to 
this parameter (usually scaling exponentially with k). 

• Range (£) - fixed the number of involved sites per correlator C'"' we 
can still choose the maximal distance among such sites for which we 
require the corresponding C'"^ to be present in the product. 

• Shape - if we are dimension higher than 1, not only the number of 
involved sites per factor, but also their shape in the lattice is relevant. 
For example, a common fashion for square lattice is work with plaquette 
correlator product [56], where every C'"^ involve the /c = 4 vertices of 
a lattice square (or a rescaled lattice square). 

At any rate, CPS are are considered successful variational tools, for their 
ridiculously small number of variational parameters, high manipulability, and 
the capability of adjusting any single correlation by adding just one suitable 
factor to the product. Choosing a suitable starting state |$°) is also a delicate 
issue for the CPS ansatz. For instance, a MPS (with modest bondlink D) 
could be a suitable choice; this lead to the construction of hybrid MPS^CPS 
Tensor Network designs, also known in literature as Renormalization Algo- 
rithm with Graph Enhancement (RAGE) 



4.8.1 CPS MPS/PEPS 

We would like to point out an interesting fact about CPS. a Correlator Prod- 
uct having a finite range i of factors which does not scale with system size 
L clearly describes a finitely-correlated state: this tells us that finite-range 
CPS can be put in relation with MPS and PEPS, or, more precisely, with 
Matrix Product Operators and Product Entangled Pair Operators (PEPO - 
operatorial version of PEPS). A result of my personal thesis work, is that 
when i is finite is always possible to represent a CPS as an analytical sub- 
class of MPS/PEPS with finite correlation bondlink D being a function of 
i. Such construction obviously depends on the involved sites per factor and 
their shape, so we are here going to present the mapping CPS — )■ MPS in 
the simplest setting: ID, and binary factors k = 2, kept at all ranges i' from 
neighbors up to an arbitrary finite bound 2 < i' < i. 
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Let us give this prescription starting from i = 1 and then increasing i: 



(4.37) 



S1---SL 



where we are assuming that C^^'-^l may depend on the pair of sites {j,j + 1} 
on which it is acting, to keep major variational freedom. Also, boundary 
conditions are small issue: they only set the upper bound L for the j over 
which the product is taken (to L for PBC, L — C for OBC). Let us draw the 



diagrammatic representation of (4.37), given by 




(4.38) 

As you see, it is possible to embed the correlator within an MPO structure 
O'-'l where the product matrices correspond to the blue boxes in (4.38), or 
equivalently 



Of" = Y.^r,s ci;^, \s){s,^A = J2 C c^f 



(4.39) 



Q,/3 



Now let us move to ^ = 2, where next- nearest neighboring correlator fac- 
tors C^^'-^'are introduced. Then the overall correlator product operator is 
expressed by: 



0[il 




(4.40) 
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It is easy to check that the circuitry of hnks and delta nodes performs the 
right index plugging at the correlator factors. Indeed, focus on C'^'-'^ in 
the picture, while it is obvious that its left link corresponds to Sj, you can 
similarly follow the right link and verify that it is exactly Sj+2- So the 
Matrix Product Operator block O^^ is also correct, as it contains all the 
inner structure (deltas) and variational information (correlators) needed. 



Equation (4.40) has a vertical pattern that is suitable to be copied and 



repeated, every time we add a layer, it is equivalent to add a new set of 
product correlators, effectively enhancing the maximal range i by one. So if 
our purpose is to describe the action of 



I L 



(4.41) 



then the corresponding Matrix Product Operator representation would read 



= E n c'S^ ) ^^Lr I n K+.^.' 1 i«i • • • • • • (4-42) 



\l'=l 



For better comprehension, let us sketch it for i = 4: 

Ob] 




(4.43) 
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The total correlation bondlink dimension used for this representation is d^, 
it should be convincing that this representation is also optimal (if no further 
restraints on the C factors are requested), since we are transferring through 
the MPO bondlink the minimal information to reproduce the factors exactly. 

In conclusion, we can establish a general entanglement bound on ID bi- 
nary Correlator Product States, which is 

5vN<^logrf, (4.44) 

characterizing the correct ID area law of entanglement if the maximal range 
£ does not scale with the total length L. 

The scheme we just presented has been originally conceived by me during 
my doctorateship work, and supported by R. Fazio and F. Becca. 

4.9 Towards hierarchical Tensor Networks 

Tensor Network architectures can be put in tight relation with numerical 
renormalization groups. In fact, consider a tensor attached to some of 
the physical indices; it can be as well interpreted as the action of a linear 
transformation acting on the local density matrix, mapping it into the virtual 
links space. As the overall effective dimension is typically reduced, a numer- 
ical renormalization is taking place. Tensor Networks entanglement bounds 
guarantee that the amount of correlation in the TN tailored variational state 
matches the entanglement that can be built via renormalization process. 

We want now to analyze and describe detailed properties of another 
renown class of Tensor Networks, corresponding to the original Wilson's 
real-space numerical RG. We are talking about Tree networks, and of course 
also their recent generalization: Multiscale Entanglement Renormalization 
Ansatz (MERA). These two network geometries share the intriguing property 
of embedding a scale-invariance in their pattern, so they are ideal candidates 
for dealing with critical models. We will classify them together as hierar- 
chical Tensor Network states because network nodes are linked according to 
hierarchical relations. 

The interest revolving around these methods, and the following in-depth 
analytical study, are such that it is appropriate to devote an entire chapter 
to describe their features and peculiarities. 



Chapter 5 

Trees and MERA 



When the Matrix Product State representation as tailored wavefunction 
paradigm version of the DMRG was reahzed and understood, it was an 
important breakthrough. But it was with the advent of Tree Tensor Net- 
works (TTN) jnU [27] and Multiscale Entanglement Renormalization Ansatz 
(MERA) [m [Tg Ha |20l 1211 [221 [28] that the computational physicists' com- 
munity started talking about Tensor Networks in general. TTN and MERA, 
similarly to MPS, have a network pattern which is not only very simple and 
highly adaptive but also self-similar. Nevertheless, the relevant difference 
between TTN/MERA and MPS is that, while self-similarity in MPS arises 
every time a new site is added, for (binary) TTN/MERA it happens every 
time the number of sites is doubled (or tripled for ternary TTN/MERA, and 
so forth). This suggests us that, as we will see, TTN/MERA bear some- 
how a resemblance to real-space renormalization processes, and also that 
they should be especially suitable for describing systems in which a scale 
invariance is emergent. Most of the results provided in this chapter are orig- 
inal work developed by me and was supported by V. Giovannetti, M. Rizzi, 
S. Montangero and R. Fazio. 



5.1 Real-Space numerical renormalization 



In section 2.1 we argumented thoroughly how a Matrix Product State can 
be explicitly built from the data of the recursive (standard) density matrix 
renormalization transformations. Let us follow a similar path now, but with 
a different starting algorithm in mind. For consistence we will work in ID. 

We are again assuming that numerical renormalization process is being 
performed, but not according to the traditional DMRG framework, where at 
every step the degree of freedom of a single site is added to the picture, and 



105 



106 



CHAPTER 5. TREES AND MERA 



the joint density matrix is then renormahzed, i.e. old block U added site — )■ 
new block {D®d D). Instead, here we are assuming that the renormalized 
degree of freedom of the old block is doubled, or equivalently, coupled with 
a copy of itself. Now we select a new state within this space (according 
to whatever benchmark we prefer), describe its density matrix, and then 
renormalize, i.e. old block U another old block — t- new block {D ® D ^ 
D). What inspires this proposal is that during the same time we perform 
a density matrix renormalization step, ID lattice sites in the real-space are 
also renormalized, in a coarse-graining fashion. 



Let us now recall equation (|2.2|), used in a standard-DMRG, while we 
! right-propagating the 
.trices = Y!j=iPj \L 
expressed through A^^^ as 



are right-propagating the scheme, and keeping D x D renormalized density 
matrices pi = J2f=iPj The basic step of the algorithm could be 



D d 

Ms I r \L 



k=l s=l 



and became the elementary tensor of MPS description. Let us write 
the corresponding recursive expansion for real- space numerical RG: if no 
translational invariance assumption is made (for instance, we could be in 
an OBC setting) we must keep track of every pf'\ where h tells us how 
many times the renormalization has been performed already, and i carries the 
information on which original sites the density matrix is actually describing. 



With these considerations, and defining p^^ = XlfLiPi i-^jl^ S^^ 



m?= Yl <M |M..)S:?®|ii/4,)g:;]. (5.2) 

fcl,fc2 = l 

where h goes from to h = logL, with L the system size, while i ranges 
from 1 to 2~^L (notice that the number of allowed i values halves at every 
layer h increase, which is exactly the coarse graining). Finally, iMj)^^ goes 
from 1 to the renormalization dimension Dh, eventually depending on h. 

You should convince yourself that with this procedure, the state is com- 
pletely and uniquely characterized by the set of A''^'^', which is a three indices 
complex tensor. Since the global state |\E'), when written in the canonical 
separable basis, has components determined by linear contraction relations 
between the A''^'^', it is a Tensor Network state. Precisely, it reads 
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where we wrote the A as Df^ x D^^i matrices. Equivalently, in the diagram- 
matic Tensor Network representation it appears hke this 




(5.4) 

as you see it is a tree graph, with branching number 6 = 2. The tensor C 
standing at the top of the structure is the only one which is topologically 
different, as it has only two links, and due to its placement is often referred 
to as hat (or root) tensor. 



From equation (5.4) the coarse graining action performed by renormal- 
izers A becomes immediate and clear: every layer of A tensors maps an 
adjacent pair of (eventually already renormalized) sites into a single renor- 
malized site, thus effectively halving the overall size of the system. Once L 
has been reduced to 2 we describe it as a simple binary state \C). 

Due to hermiticity of every p^^^ we would like for their respective eigenvec- 
tors iMj)^^ to be an orthonormal set, at every h and i. As for MPS, where a 
similar requirement lead to a gauge symmetry breaking, this restraint trans- 
lates into a condition that every A must satisfy, namely 



6j, 



J2 



E 

ki,k2 



W{h, 



(5.5) 



where * stands for complex conjugation. In other words, every A, read as a 
Dl^_-^ X Dh matrix, must be (left-) isometric, i.e. A"''A = 1. This is indeed 
a gauge symmetry breaking. Truly, this is exactly the peripheral gauge we 



defined in section 4.5.1, when the nucleus 11 corresponds to the hat tensor 
C. We remarked that, for a given Tensor Network with no closed loops (as 
a tree graph is), it is always possible to find the gauge transformation that 
maps it into the peripheral gauge, no matter the starting state: this tells us 



that the isometricity condition (5.5) for A carries no loss of generality at all. 
You can guess that the binary character 6 = 2 of the tree network in 



(5.4) is due to the fact that we renormalized just two copies of the old block 



density matrices into a new one. Of course, mapping an arbitrary number 
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b of copies of the old block into a single one leads to tree Tensor Networks 
with the corresponding branching number b. For example, we could have a 
ternary tree Tensor Network when 6 = 3: 




(5.6) 

where the total number of layers h would he h = \og^ L. Throughout this 
chapter we will present various results, deriving our calculations in accor- 
dance to the binary tree Tensor Network (2TTN) design. We want to remark, 
however, that most of the claims and statements extend naturally to TTN 
with higher branching number: typically one has just to substitute 2 with b 
where appropriate (usually as a logarithm basis) to get the right result. 

Moreover, notice that we disregarded any assumption on boundary con- 
ditions so far. Indeed, unlikely from Matrix Product States, where different 
boundary conditions lead to substantially different network topologies (no 
closed loop in OBC, but with closed loop in PBC), in TTN representations 



of a state, like (5.4) and (5.6), the network topology is insensitive to the 



presence of a boundary. This, as we will see, will naturally bring a unique 
and well-built definition of the thermodynamical limit. 



5.2 Tree network entanglement instability and 
the introduction of MERA 

If we wish for our real-space self-similar Tree network to be actually capable 
of describing ID critical states, we must fist check that the amount of entan- 
glement it can hold satisfies the typical entanglement area-law violation: 

SYNipe) = -\ogi + const. (5.7) 

where pe is the density matrix of i adjacent sites, and in PBC conditions. 
Also we are assuming that we are approaching the thermodynamical limit, 
or at least that the region i is too small to capture finite size effects: I <^ L. 

In order to evaluate upper bounds to the entanglement of a TTN, we are 
going to employ the Tensor Network entanglement arguments we provided 
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in section |4.2[ There we found that the entanglement of a partition is 
bounded by ~ where x is the minimal number of broken links needed to 
separate the network graph into two subgraphs, respectively containing Ti 
and its complementary F^. Let us perform the count for binary TTN. First, 
notice that, since the Tree network design is not translationally- invariant de- 
fined, we will expect that 5vn(P{^i/i+^}) will not depend only on the number 
of sites i in the interval, but also on its placement ii. 

Now we proceed recursively with layers. If ii is odd, then I have to cut no 
link in the lowest layer, if it is even, i cut one link; either way, i move up one 
layer and site ii is mapped into renormalized site [^i/2] (where \x] is the 
smallest integer which is larger than x), and we can repeat the procedure. 
The same argument holds for the other region boundary ^2 = + ^, even 
though the link breaking is needed when £2 is odd and not when it is even. 
In conclusion we have that iSvN(p{fi/i+£}) < #cuts ■ logD, where 



1 < #cuts <2h = 2 log,, 



(5.8) 



depending on £1. For clarity, let us show an example of the described proce- 
dure when L > 16 {h > A), are we are considering P3,ii: 




(5.9) 

where we have drawn only a branch of the whole Tree network. In this exam- 
ple the entanglement bound is given by ~ 4 log D (actually by ~ [3 log D + 
log d] since one of the links to be cut is a physical one) . 



Result (5.8) is quite relevant: it tells us that it is true that in a TTN 
state one can choose a sequence of subsystems of growing size £ exhibiting 
the correct area-law violation scaling, but it is also true that it exists another 
sequence which is nearly separable. In other words, we can state that the 
entangling capacity of a TTN fluctuates widely with translations, ranging be- 
tween almost no entanglement and critical entanglement. Such translational 
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instability is embedded within the nature of Tree networks themselves, and 
must be handled with care. In particular, in order to work around this en- 
tanglement instability, and representing a variational quantum state showing 
a more smooth area law violation, two methods are adopted: 

Incoherent translational mixture - This path focuses on calculating 
the translational average every time we wish to acquire an expectation value 
on the TTN state. In practice, let 6{f^}^ be an observable having support 
on sites then what we are actually interested in is 

1 ^ 

<®) = zE(%+mJ- (5-10) 

By doing so we always integrate out translational fluctuations, as we are 
considering solely the zero Fourier mode. This operatorial average procedure 
is theoretically equivalent to consider a state which is the incoherent mixture 
of all possible translations of the original TTN state. Namely we define 



1 ^ 

j^P{t+t^}^- (5.11) 



for any choice of the support {£a}a and we have represented a state, which 
is translational by construction: P{^^}„ = P{e+ia}a- Notice that we are not 
considering the coherent superposition of the TTN state translations, because 
the interference graphs can not be contracted efficiently due to the presence 
of several closed loops. We can now combine the entanglement argument of 
Tensor Networks with the concavity property of Von Neumann entropy, and 



obtain an overall bound on the entanglement for state (5.11): 



5vn(p^) <log2MogD (5.12) 

which satisfies the right logarithmic behavior, and is perfectly smooth under 
translations. 

Adding new tensor elements: MERA - The idea of MERA arose 
from the request of saving some information on correlations shared by neigh- 
boring sites scheduled to be renormalized into different, separate blocks. For 



instance, in (5.4) site 7 and site 8, despite being neighbors, are not going 
to be renormalized together until the hat is reached: this means that most 
of the entanglement they share is likely to be lost in the renormalization 
process, and poorly described by the corresponding Tensor Network design. 
Therefore, the scheme proposed by G. Vidal [TTj consists into performing a 
(unitary) operation X coupling these sites, whose purpose is to disentangle 



5.2. TREE ENTANGLEMENT AND MERA 



111 



the two as much as possible before renormahzing them separately. This way, 
the information concerning the original entanglement of the pair is stored 
within this operation X, and renormalization needs not to concern about 
it. Disentangling operations and real-space renormalizators are then applied 
alternately. 

It is clear that, since unitary disentanglers are linear operations, quantum 
states achievable by this process are again Tensor Network states. They are 
known as Multiscale Entanglement Renormalization Ansatz states, and show 
a hierarchical and real-space self-similar pattern just like TTN states. Here 
we present a partial frame of a binary MERA Tensor Network 




(5.13) 

where for every layer h and horizontal position i, tensors are chosen to be 
unitary/isometric, i.e. X'^X = 1£)2xd2 and A''"A = Idxd- Of course, the 
corresponding MERA version can be modeled on a Tree Network of arbitrary 
branching number; in literature, both binary and ternary MERA have been 
considerably used for simulation purposes. For completeness, let us draw the 
diagrammatic pattern of a ternary MERA: 




(5.14) 

We will show now that introducing the new Tensor elements X plays a fun- 
damental role in regularizing the entanglement scaling law under translation. 
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even though MERA are again, hke TTN, a non-translational network design. 

Let us assume, as before, that we are to evaluate the entanglement en- 
tropy of estimating an upper bound via the Tensor Network link-cut 
argument, and proceed recursively layer- by-layer (a MERA full layer is com- 
posed by stacking together disentanglers X'^''^ and coarse-grainers 

AIM). 

For every layer, and each one of the two region boundaries {ii and i + ii), 
I have to cut either one or two links, depending on the region location with 
respect to the MERA network geometry, as you see from this example 




(5.15) 



This leads to a bounding function for the entanglement i5vn(p^i/i+^) < 
#cuts ■ logD, where 

2 logb i < #cuts < 4 logb i, (5.16) 



depending on ii. It is clear from (5.16) that entanglement-bound fluctuations 



are present even in the MERA case. However, differently from (5.8) they are 
hardly a problem: indeed even in the worst case scenario iSvn(p^i,£i+£) < 
21og^£ ■ logD: the entanglement is ruled by a logarithmic violation of the 
ID area-law, as we wished. This is the main reason why MERA are kept in 
high regard in ID critical systems simulations. 

Despite the enhanced accuracy given by disentangling elements of a MERA, 
we would like to show, throughout this chapter, that MERA and TTN rep- 
resentations manifest really a common behavior, as they can both easily cap- 
ture critical properties of strongly-correlated states, e.g. in terms of critical 
exponents, primary fields, and so forth. Acknowledging that, intuition sug- 
gests that renormalizer operations A are those responsible for keeping track 
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of long-range properties of the TTN/MERA state, while disentanglers X are 
mostly used to adjust the local variational structure. To avoid translational 
invariance-breaking issues, we will make large use of translational averages of 
operators and correlators in the following sections, thus substantially adopt- 



ing the incoherent mixture paradigm of (5.12) whenever possible, even for 
MERA. 

Similarly, notice that TTN can be definitely seen as a subclass of MERA. 
They actually correspond to MERA networks where all the disentangling 
tensors are set equal to the identity, i.e. X^^'^^ = 1£)2xd2 V/i, £. 

One additional remark: we stated that Tree networks are suitable for any 
boundary condition we might inquire, this is due to the presence of dense 
graph frontiers in their geometry. On the other hand the MERA Tensor 
Network has a natural attitude for preferring Periodic boundary conditions. 
Indeed, as you see from (5.13) and (5.14), each site has a network path 



through the lowest layer linking it to both its neighbors, recalling somehow 
the valence bond picture idea for Periodic MPS. Precisely, one can represent 
the whole ID binary MERA state as follows: 




(5.17) 

(figure taken from ref. [23]) where the ring of physical sites stands at the 
edge of the circle. 3-link nodes are A tensors, 4-link ones are X tensors, and 
the tensor in the very center of the diagram is a 4-legs hat tensor IC4). You 
can easily check that (5.13) is exactly the local pattern of (5.17). Of course. 



a generalization of MERA for OBC can clearly be designed by adjusting its 



geometry, we will treat this point in section 5.10 
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5.3 Causal cone property 

A remarkable feature, important both from analytical and computational 
point of view, shared by TTN and MERA alike is the so-called causal cone 
property. It is a generalization of the peripheral gauge paradigm, allowing 
us to contract at no expense whole branches of the Tensor Network, that 
extends also to MERA geometry, despite having closed loops in their graph, 
thanks to the isometricity requirements for X and A tensors. 

Assume we want to achieve the expectation value over a ID TTN/MERA 
state for an observable Q{ii..£2} whose support is an interval of sites, say 
{^i..^2}- Then, by adopting the standard network contraction method, one 
performs (*|9|^') = E{r},{s} '^r}(^^i • --ri^lQlse^ . . .s^j)?^^}- But now, con- 
sider the lower layer of disentangler X^^'^' (for a MERA), which are unitary 
by requirement X'^X = 1. Those X whose (both) physical indices do not 
belong to the support interval {ii..i2} are contracted with their adjoint and 
automatically vanish. Precisely, we got rid of a whole layer of disentanglers 
except for those directly connected to the support. The same argument holds 
for renormalizers A belonging to the lower layer, which disappear alike, as 
sketched in this diagram: 




(5.18) 

As we eliminated an entire (double) layer from the Tensor Network state and 
its conjugate, we could say that we mapped the original expectation value 
problem into a new one: 

(*|e|^) = (^'|^|,,..,,}(0)|^') (5.19) 



where the quantum state |\E'') is the original TTN/MERA state without the 
bottom layer, and the new observable Aii-^^.j^} (O) is the big composite pink 
tensor of equation (5.18 b). Acquiring ^{^^..^j} (©) starting from 9, X and 
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A tensors is a finite operation, whose complexity scales with the size of the 
support £2 — ^1, not with the size of the system L. Moreover, when mapping 



9 — )• A{e-^,,e2} (©) the actual support typically shrinks: you see how in (5.18) 
we map a 5-site observable into a 4-site one. We can obviously refrain this 
argument, every time mapping up the effective observable in a recursive 
fashion 

ew^0[/^+i] = ^w^^^,^^(eW) (5.20) 

until we reach the hat tensor \C), for which it holds (\I^|6|\I^) = {C\Q^^~'^^\C). 



Each mapping contraction (5.20) has a limited computational cost, and we 
have to perform it a number of times equal to the total number of layers 
h, but we know that 2^ = L, the total system size. In conclusion, for bi- 
nary TTN/MERA the computational cost for evaluating compact-support 
observables scales as 

#cost oc log2i^. (5.21) 

Obviously, tensors canceled out by isometricity relations can not influence 
in any way the result of the expectation value, so all the physical properties of 
sites in {£1 ... £2}, and so their reduced density matrix, can be determined by 
a small subset of tensors in the network. Namely, graph nodes we can reach, 
starting from physical links {£i...i2}, by means of sole vertical network 
propagation (only moves that increase the layer index are accepted). If we 
want to use a language formalism borrowed from relativity, those 'influent' 
tensors form the causal cone of sites {£1 . . . £2}- Here we show you an example 
for a binary MERA with L = 16 sites: 




the yellow region is the causal cone of sites {4, 5, 6}, it contains only 4 disen- 
tanglers X, and 5 renormalizators A, out of 25 tensors in the whole network. 

Causal cone 'horizontal' sizes (width) are dynamic while moving across 
the layers, as already pointed out. In practice, given support size and place- 
ment of the pre-mapping observable, by simple diagram contraction rules one 
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can determine the size of the mapped observable, which depend also on the 
chosen network geometry and branching number b. Precisely for 6 = 2: 



binary TTN 



2i-l 
2£ 



i if ii odd 
+ 1 if £i even 



binary MERA 



2^ - 1 ^ + 1 

^+1 if^iodd (5.23) 
? + 2 if £i even 



2i 



We immediately see that for binary TTN geometry, a cone width of 1 is 
stable through the causal path. Accordingly, a size of 3 or greater shrinks in 
an almost-exponential fashion with layers. Finally, an observable supporting 
2 adjacent sites may either be mapped into a one-site observable, or a two- 
site one, depending on how these sites match the network structure: we will 
say, then, that 2-sites is a metastable size, as it may collapse to 1 but not 
necessarily, nor immediately. This remains true for higher branching number 
b TTN states. 

For a MERA state it can even happen that we have to grow the causal 
cone width, e.g. if the observable was one-site is mapped into a two-site 
one. In binary MERA geometry, 3 is the only stable size, while 2 and 4 are 
metastable. When moving to a ternary MERA geometry (or a higher branch- 
ing b one), we see that these characteristic widths are smaller: precisely size 
2 is stable, sizes 1 and 3 metastable (but 3 becomes unstable if 6 > 3). 



5.4 Ascending and descending maps 

In this section we will focus on those maps that propagate operators upwards 



trough the layers, like the one of equation (5.20). According to their nature, 
we will call them ascending maps, and relate them to completely positive 
maps described in appendix [A} 

For simplicity, let us start with a one-site operator 6^ acting on a binary 
tree state. We acknowledged that size 1 is a stable causal cone width for this 
TTN geometry, so the mapped operator A{i}{Q) will still be one-site. Let 
us express analytically the action of this mapping, which in turn depends on 
the parity of i 

^ ^ ^ f ^L(e) = At(e®i)A iff odd ^ ^ 

Asn{Q) = I ^ ^ / ^ 5.24 
^ ^ [ ^i?(e) = At(l®e)A if £ even 
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where the subscript L/R refers to which of the two lower hnks of the node 
tensor A (left or right) is touching site i. Of course, tensor A also depends 
on the layer h and horizontal position i that locate its placement within the 
network, and thus similarly will do maps Al and Ar (however we shall treat 
implicitly this dependence to avoid carrying on too many indices). 



It becomes immediately clear from (5.24) that both Ar and Al maps 
are Completely Positive and Unital. This is due to the fact that we already 
expressed the Kraus expansion of the mapping, indeed A* ^ form a Kraus set 

over k of D X D matrices A^^^ = j fcK) 01' ^^"^ another set over j of 

matrices A^^' = J2i k^] k\'''){^\- Therefore 

^^(e) = ^Af]^eA[^] and ^,^(0) = Af ^GAf (5.25) 

k j 

but due to isometricity condition on renormalizators, we have 

^ Afl^t^J = ^ Af Ufl = AtA = ^^(1) = AlW = 1, (5.26) 



which is exactly the Kraus set requirement (A.l). Interestingly enough, it 
is relevant that this map should be completely positive; the reason becomes 
clear one we introduce the formalism of descending maps for density matrices. 
In fact, since 9 has 1-site support, the reduced density matrix of that site is 
the sole responsible for determining expectation values: 

(^|e|^) = Tr[e-p,]. (5.27) 

The causal cone argument tells us that this is also equivalent to contracting 
A{£}{Q) with the TTN state without the bottom layer, on which A{e}{Q) 
acts again as a single site operator, so that 

{^'\A{emm = Tr [A{em ■ p',] = Tr [0 . 1?m(p^,)] , (5-28) 

where p'^, is the reduced density matrix of |\I'') at site i' = and V^iy is 

the map adjoint to A{e} with respect to the trace scalar product for matrices 



{A\B) = Tt[A^B]. It is clear that since (|5.28l) holds for the whole algebra of 



1-site observables G it must necessarily be that 

p, = V{eM,). (5.29) 



Luckily T>{£j, being adjoint of a CP-unital map, is a completely-positive trace- 
preserving (CPT) map, which takes density matrices to density matrices and 
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thus (5.29) is perfectly meaningful. As pointed out in ref. [22], applying V^i-^ 



can be actually seen as the action of a quantum channel pQ. The equiva- 



lence in (5.28) can be interpreted as some TTN version of the Heisenberg 

we can map ascending-wise operators, or 
the expectation value is achieved either 



f-)- Shrodinger scheme duality: 
descending-wise density matrices, 
way and the number of contractions is the same. Similarly to A{e}, T^{e} also 
depends on the parity of target site £, namely 



I)^(p')=Tr2 [Ap'At] =5^Af]p'A 



©^(pO^Tri [Ap'At] =5^Af p'aJ- 



[R] J \ [R]^ 



(5.30) 



both CPT, and expressed in Kraus formalism. Performing the algebraic 
contraction of any of the maps defined in (5.24) and (5.30) requires a number 
of elementary operations scaling like ~ D^, as you can easily guess from the 
following diagram representation: 



A 



R 




V 



R 



(5.31) 

We presented the ascending/ descending map formulation for one-site op- 
erator/density matrix, but we are not surprised that it extends naturally for 



any size, according to the causal cone width rules (5.23). Clearly, ascending 



maps will typically reduce the size of the effective operator, and at the same 
time, descending maps with a given target size will be function of density 
matrices with smaller size. Let us consider the case of two adjacent sites 
+ 1}, depending whether ii is even or odd we have: 



P{2i-i,2e} 

P{2l,2l+l} 



5(p^) = Ap^At 

^^i?®^^L(P{£,f+i}) 



(5.32) 



where S is the CPT map obtained by just multiplying the density matrix by 
A and K^] and notice that it holds Vl{-) = Ttr[S{-)] and Vr{-) = Ttl[S{-)]. 



The map appearing in the bottom line of (5.32) is separable, so it can only 
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decrease entanglement, and likely it will since it is not an invertible trans- 



formation and typically satisfy the mixing property (see appendix A.3). Let 



us sketch the diagrams for the descending operators of (5.32) as well: 



s — 




m i 









L 




(5.33) 



the ascending ones are simply defined by taking the adjoint of (5.33). Due 
to the fact that TTN has a renormalization flow which is quasi-local at every 
step, it is easy to check that maps S, Vl and are the only ingredients we 
need to build all the maps for any given size. For instance, let us write the 
generic descending map for the binary TTN, in order to obtain p|^^ 



odd length 



P{£i,ii+2£} — 5 



' ^ ® 1^L{p'{iej2],iei/2]+e}] 

{£i/2,£i/2+£} ) 



for £i odd 
for £i even 



even length 



P{hA+2£-i} = S® ... ® S{p[^f^/2],\e,/2]+e}) for ii odd 
® 5 ® ^^ij(P{^,/2/i/2+f+i}) for 



Ci even 
(5.34) 

We agree that working with all these even versus odd possibilities looks defi- 
nitely messy and confusing. To work around this issue (even though partially) 
we suggest to employ the translational average incoherent translational 
mixture formalism. That is what we will discuss in the next section. 

Ascending/descending maps for MERA - Before moving on, it 
should be pointed out that, as we developed an in-depth formalism of com- 
pletely positive ascending and descending maps for binary tree network, this 
can be similarly done for other TTN geometries and for MERA as well. 
Precisely, in MERA geometries, maps also embed the local action of dis- 
entanglers X as well as renormalizers A. As an example, let us write the 
descending maps for 3-site density matrices (recall that 3 is the stable causal 
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cone width) V^f-^ and V^f-^, respectively reading 





defined in this fashion so that 



/2A/2+2}. 



„MERA 
/^{£i/i+2} 



'T-.MERA^ „MERA 

V/^MEB 

'3-s>3,fll/^{L4J,UiJ+2}/ 



V 



MERA/ MERA 



for £i even 



for i-i odd 



(5.35) 



(5.36) 



Extending the descending map formahsm to density matrices of any size is 



performed similarly to (5.34) for MERA, but with a main difference: the 



maps can not be any longer expressed as tensor products of local objects due 
to the presence of disent anglers. 



5.5 Horizontal homogeneity and 
translationally averaged maps 

In TTN/MERA Tensor Network, interesting new properties arise when we 
adopt the additional assumption that tensors belonging to the same layer 
h of the network structure (and of the same type X, A for a MERA) are 
identical, that they are actual copies of the same variational tensor, appearing 
times at different nodes in the network. We will call this requirement 
horizontal homogeneity, because we are still letting the tensors belonging to 
different layers free to be variationally different. 



In section |3.2| we showed how tensor homogeneity in periodic Matrix 
Product States is a condition strictly related with translational invariance 
symmetry. Neither for Tree networks nor for MERA, the relationship be- 
tween homogeneity and translational invariance is so sharp and well-defined, 
since the network geometry itself is visibly not translational. Nevertheless, 
horizontal homogeneity is not only the best approximation for emulating a 
translational invariance, but also, and more importantly, allows us to build 
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a well-defined ascending/ descending map formalism for the translational av- 
erages -H- incoherent translational mixture scheme. This will let us write 
simplified versions of the (say, descending) map definitions we provided in 



last section, namely (5.30) (5.32) and (5.34), as they will depend no longer 



from the placement index ii, but solely on the size £. 

Then let be the one-site density matrix of the TTN, averaged over 
translations 



Pi 



2M 



+^-1}' 



(5.37) 



ii=i 



where now the layer index /i is counted starting from the hat and moving 
down (i.e. n = h — h), for comfort. This will prove a labeling advantage 
when we will discuss the thermodynamical limit n — >■ -|-oo. 

As before, we start from smaller sizes and then grow with increasing i, to 
enhance immediateness and comprehension. Consider then , by exploiting 



(5.30) we can write 



^1=1 



P^ ~ Ii' 2^ ^{^i} ~ 2^ ^ ^{2^-1} ^ ^{2^1} 
- - ^1 = 1 

= ^E^(^^'H''1m')+^"(''S«)) (5-* 



where and are those descending maps defined in (5.31 ) when tensor 



A belonging to layer p is being used. Then by linearity of the descending maps 
we can re-sum the argument density matrices to obtain again the appropriate 
translational mixture 



Pi 




2''-i 

— V 0^''- 

L-\ P{IX 



} 



Pi 



V^^Hp\-''^\, (5.39) 



where we defined the average map T> = (T>l + T>pi) /2. This new map is again 
CPT, as it is the equally- weighted mixture of the two original maps, its 
Kraus decomposition is given by the union of the two Kraus decompositions 



of T>L and Vp^ (every matrix multiplied by a 2 factor). Eq. (5.39) yields 



to a single layer-to-layer recursive relation 



instead of the 



whole family of equations (5.29 ) and (5.30 ), thanks to horizontal homogeneity 
requirement. 
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Let us move to two adjacent sites, a causal cone size which is metastable 
for TTN geometry: its stabihty will be regularized when performing trans- 
lational averages, as we will see. Then 



P2 



2^-1 



(5.40) 



Watching this result, one can argue that in order to achieve pg^' by descending 



recursive relations, we should need both ^' and pg'^ ^' 
contains by itself all the information concerning Pi^~^': 



definition (5.37), which tells us that Tri[p2 



But indeed, 
this is obvious by 



present consideration allows us to write 



T^'2X2{pr'') where 



The 



n 



[m] / 

2^2 1 



2 ^ 



+ (Tr,[.]) + (Tr^[.]) . (5.41) 



cos^ 6 and sin^ 6 are just two free real parameters which must be both positive, 
to guarantee complete positivity of ^^2^.2 5 ^.nd whose sum must be 1 to ensure 



trace preserving property. Then for any real angle 9, eq. (5.41) performs the 
right mapping. 

By natural extension, we can generalize descending mappings for any 
density matrix size, which can be always written in terms of iS, T>l and Vr. 
We report here the result for completeness: 



T)i^2i-i = ]^{Vr®S® ...®S + S 



S®Vj 



(5.42) 



where Tr^ means tracing away the rightmost site times cos^ 9 plus tracing 
out the leftmost site times sin^ ^ (angle 9 arbitrary), as in (5.41). All these 



maps are convex combinations of the maps previously defined, so they are 
indeed CPT. 

Ascending maps can be naturally defined in a translational invariant fash- 



ion as well; it is easy to see that the composition relations (5.42) hold for the 
adjoint maps A2i-2^t and A2t-2^t identically. 
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5.6 Complete homogeneity and 

definition of thermodynamical limit 

We mentioned that the main purpose of introducing Tree and MERA Tensor 
Networks is to simulate strongly-correlated systems. It is clear, though, that 
this ansatz would be more suitable the more we are capable of reproducing a 
conformal symmetry, where translational invariance, and more importantly 
scale invariance hold. 

In previous section we dealt with translational invariance by merging 
together the horizontal homogeneity requirement and the incoherent trans- 
lational mixture framework. Now we are going to force scale invariance into 
the system by assuming complete homogeneity: i.e. requesting that ev- 
ery renormalization tensor A in the network is identical to the others, even 
throughout the layers. We will see that this assumption will naturally lead to 
the definition of a thermodynamical limit, in which every physical quantity 
is perfectly controlled thanks to isometricity relations, and where the scale 
invariance becomes manifest. 

This can be explained heuristically, by considering the family of tree 
graphs. Every tree graph has a self-similar pattern, but only for an infinite 
tree the self-similarity becomes exact: every branch is topologically identical 
to each of its sub-branches because they are all infinitely-long. Such sym- 
metry, though, would be broken if one were able to distinguish somehow the 
graph nodes, thus the assumption of complete tensor homogeneity. 

Let us then consider the sequence of completely homogeneous TTN states 
I^M^A)), defined solely by a single tensor A, repeated at every network node, 
and indexed by the total amount of layers fi. We are going to characterize 



the limit of this sequence for /i — > cxd. As we discussed in section (3.4) the 
thermodynamical limit state is defined by the family of density matrices 
for every size £, with the requirement that elements of this set undergo the 
correct partial trace relation. In our TTN setting we are first considering 
translational averages, and then approaching the limit, so that 

= lim pfl (5.43) 
We will again proceed starting from size 1. Horizontal homogeneity tells us 



that (5.39) holds, but at the same time, vertical homogeneity implies that V 



does not depend on the layer we are considering. Therefore 

= V (p[^-^l) = (pj^-'l) = ... = Vf^ (pS^^^*!) , (5.44) 

where pf is the translationally averaged one-site density matrix of the hat 
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tensor. Then we have that 



= hm 

fl—^OO 



Pi ) 



(5.45) 



i.e. one-site physical properties of the thermodynamical hmit state are given 
by applying V to the hat tensor. We will now make an additional assump- 
tion, which is indeed a constraint on A but formulated it in terms of V: we 
will require that D is a mixing CPT map (see A.3). This condition on the 



eigenspace decomposition of V not only states that V has a unique fixed 
point pf, but also that pj is the only attraction pole, so that every state is 
mapped into pf after infinite applications of V. Then is exactly the fixed 
point of V: 



v{pr 



(5.46) 



As you guessed, normalization is kept correctly under control at every finite p, 
and thus in the limit of an infinite amount of layers p. Indeed, we argumented 
that isometricity of A is equivalent to the peripheral gauge, for which the state 
normalization condition is automatically dumped onto the hat tensor alone, 
1 = (v['M(A)|\E'['^l(A)) = (C|C); trace preservation property of descending 
maps does the rest. 

Let us proceed to size 2: the two- adjacent sites density matrix in the 



thermodynamical limit p^ should satisfy (5.40) for p — )■ oo, which reads 



pT = lsipr 



1 



Vr ® Vl {P: 



oo > 
2 J 



(5.47) 



This is a quasi-recursive relation which has also a direct dependence upon 



p5". In particular by nesting (5.47) within itself, we can rewrite the recursive 
expression into a series, namely 



pT 



,T=0 



Vt 



(5.4J 



where the expression inside the square parentheses is a CPT map, because 
the positive weights 1/2^+"^ sum to 1. At the same time, we want our pf^ to 



be consistent with equation (5.41), which becomes a fixed point equation in 
the thermodynamical limit: 



pT 



T^2^2 {P' 



oo ^ 
2 J 



1 



+ ^5(Tri[p-]) + ^5(Tr2[p^]). (5.49) 
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In order to show that (5.47) and (5.49) are actually compatible equations 
(for every 9), it suffices to prove that Tri[p^] = Tr2[p^] = ■ To show this, 
simply consider the partial trace, say with respect to the rightmost site, of 
(5.47). Then we get 



Tr2 [p; 



(5.50) 



but Vl{p^) 
equal to 2p'j^ 



[2V — Vji\{p^) and due to the fixed point property is also 
T^r{pT)- Therefore we can write 



Vn{pT-^^2 [pr]) = 2(pr-Tr2 [p^]) 



(5.51) 



which is an eigenvalue equation. But since a CPT map has spectral radius 
1, it can not have 2 as eigenvalue; thus the only solution of (5.51) is pf^ = 
Tr2 [p2°]- Similarly, we can trace out the left site in (5.47) and obtain the 
other equality: pf^ = Tri [pf]. This guarantees that p^ is unique and well- 
defined. 

Thermodynamical limit density matrices for sizes 1 and 2 are the only 
ones that require dealing with a fixed point equation to be achieved: all the 
other p^3 can be directly calculated by applying a finite number of maps to 
p^. Precisely we are referring to size- increasing descending maps of equation 
(5.42), that not only generate the whole family of pf , but ensure that this 
family satisfies the partial trace requirement. For example: 



3^5 



(5.52) 



and by partially tracing this equation, say on the two rightmost sites, we 
obtain 



Tr4,5 [pT] = ^^2^3 (Trs [V^^sipT)]) = 

I?2^3 o V2^2 {p'^) = V2^3 (pT) = pT, (5-53) 

which is the right consistency check. 

In conclusion, we characterized properly and completely the thermody- 
namical limit TTN, whose uniqueness is ensured by the mixing requirement 
of maps T>i^i = V and P2-5>2- Also notice that in this limit, any residual 
dependence on the hat tensor \C) that may linger at finite sizes, vanishes. 
We can interpret this consideration as a hint that A is the only responsible 
for capturing, and keeping track, of all the bulk properties of |\E'(A)). We 
will further argument this claim in the next section, where we will calculate 
two-point correlators and show that the manifest critical behavior of TTN 
states is ruled by spectral properties of maps, and thus by A. 
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5.7 Correlations and criticality 



A critical ground state is a scale-invariant state, whose unmistakable signa- 
ture is a power-law decay rate, with the distance, of two-point correlation 
functions. We will now investigate such correlations within a TTN state, to 
prove its criticality. The translational framework scheme will allow us define 
correlations depending on the two-point distance alone and not on the loca- 
tion; and then we will drive the results towards the thermodynamical limit, 
to better match conformal symmetry. 

Then, we start defining a correlation function similarly to what we did for 



thermodynamical MPS in section 3.5, but now we average over translations 



(5.54) 



for any pair of single-site observables 0,0'. We start by probing this corre- 
lator on the finite but fully homogeneous TTN state |^[^1(A)), with H being 
the total tree graph height. Then it is clear that 



^{'^1(0,0') =Tr (0®0') ■ [a 



(5.55) 



which extends by linearity the definition of two-point correlator for any two- 
site observable F, not only tensor product ones. The pair of two-sites density 



matrices in expression (5.55) is defined as: 



£o=l 



1 

2^ 



2M 



5Z '^^{'^0 



+l...£o+^-l} 



2M 



2^* 



4=1 

io} ^ P{io+i} 



(5.56) 



Notice the difference between these two objects: a\^^ is the reduced density 
matrix of two sites at distance i, averaged over translations. ri^\ instead, is 
the translational average of the product of one-site density matrices, whose 
sites stand exactly at distance i. Honestly, ?7^^' is a separable density matrix 



by definition (|5.56|), but not necessarily a tensor product matrix, i.e. i]^'^ ^ 



One can say that r^t^l keeps track of classical correlations between 
sites, but it is free of quantum entanglement. At the same time we have that 
Tri[af]] = Tr2[ai^^] = Tr^r^f'] = TT2[vf^] = p^^ for any t 

By exploiting the descending translational map formalism of eq. (5.42), 
we can as well develop recursive relations for the quantities ct|'^^ and ri'i^ 



For 
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simplicity and clarity we will write only the equations for an even distance 
as follows: 



a. 



21 



-0 



a 



and ^^}=■0[^\ 



where the 2-sites descending CPT map '0 is defined as 



1 



-0= - {Vl®Vi^^Vr® Vr) 



(5.57) 



(5.58) 



Notice that the propagation map for the two density matrix elements is 
the same '0, although it applies to different arguments, so the difference is 
nonzero in general. In fact, assume that we are calculating the averaged 



correlator (5.54) at a two-point distance £ which is a power of 2, say 2^ 



cg2.(r) = Tr [r faj^-^i - yir"^ 



(5.59) 



where, by definition (5.56), we have (t\ = P2 . On the other hand, we have 



that Tfi satisfies a recursive equation formally similar to (5.40), but with a 
different inhomogeneous term 



(5.60) 



where for r^Q^ it holds rj^^ = 0{ri^ 



After setting up all these ingredients, we can drive equation (5.59 ) towards 



the thermodynamical state easily, where every element gains a well-defined 
limit. Precisely: 

4Tkr) = Tr [r-;E^'^(Aa)] , (5.6I) 

where Aa is a null trace matrix, given by the series 



00 ^ 



,T=0 



{S{p'^)-VL®VR{jr)). 



(5.62) 



Here we assumed that the CPT map ^satisfies the mixing requirement, with 
jj' = lim^_>.oo ?7cr^ being its unique fixed point; clearly if is mixing, then V 
is mixing as well, since Tr2[^(-)] = ^^(■)- 



It is now worth to analyze result (5.61 ) and derive some interesting conclu- 



sions involving this correlator. Immediately we notice that when increasing 
the distance to infinity the correlation function drops to zero as it should. 
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since infinitely distant sites must not be entangled. This can be checked by 



accepting that ^is mixing and using relation (A.5); then we get 



lim 4?kr) = lim Tr \T ■ 0'^ (Aa)] = Tr [T Tr[A(T] = 0, (5.63) 



because Act is a traceless matrix by definition. 

To investigate how the correlator scales at finite distances, we first switch 



to the ascending map formalism for (5.61); this is done by just taking the 
adjoint map ^ of the superoperator Then 

c'Tkn = Tr [X^r) ■ Aa] , (5.64) 

where ^ is completely positive and unital. The linear map is not nec- 
essarily diagonalizable, nevertheless it exists at least one eigenoperator 
for each of its eigenvalues A^. Then when evaluating the correlator on that 
eigenoperator we have 

4T^(r„) = Tr {A^ (r„) • Aa] = Tr [r„ ■ Aa] = A^ '^TK^o). (5.65) 

Now recall that q is the logarithm of the distance, which leads to 

l(r„) = A^°S2^ '^f" = i^- (5.66) 

We obtained an exact power-law decay for out fixed distance two-point corre- 
lator, with the critical exponent = log2 Aq,. This is one of the main results 
I achieved during my doctorateship study an infinite homogeneous tree 
Tensor Network defines always a critical state, where the two-point criti- 
cal exponents are the logarithms of the spectrum of and their respective 
eigenoperators correspond to primary fields, i.e. exactly-scaling operators. 
Notice that since |Aq,| < 1 as CPT maps are contractive, the critical 



exponents have always a negative real part 3?(^q,) < 0. Thus, although (5.66) 
might oscillate, it always decays in modulus, never explodes. 

For a generic observable F, we can write a formal expression of its cor- 
relation function which exploits its expansion in the generalized eigenbasis 
{a, d, w} of 0, like we did for ( |3.17[ ). This allows us to write 

£f 1(F) = ■Y.'^t'''"^i^og,i) (5.67) 

a d,w 

Where V^^^ are polynomials of degree x, which eventually arise from the Jor- 
dan block structure of 0. As you see, the power-law always dominates the 
logarithmic part and rules the physical critical behavior of the thermody- 
namical TTN state at long distances. 
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a 



Figure 5.1: Critical exponents ^q, (in modulus) for an infinite Ising chain, 
calculated with homogeneous MERA algorithm. Here a bondlink D = 4 was 
used. The red dashed lines report the theoretical values. The graphic was 
kindly contributed by M. Rizzi and S. Montangero |29j . 



Criticality of MERA - The results and observations we just presented 
hold for a MERA geometry as well, although since size 1 is not a stable 
causal cone width, it is necessary to involve bigger-sized observables and 
density matrices. Precisely, in a binary MERA, even if we started from a 
1+1 (two-point) operator, after applying few ascending maps, we increase 
its support size, typically mapping it into a 3+3 operator (three adjacent 
sites, and other three adjacent sites, these two parties standing at arbitrary 
distance), and after that size becomes stable. According to this framework, 
critical exponents of 3+3 correlators (and thus also of 1+1 correlators, which 
are a subclass of the 3+3 ones) are given by the logarithms of the spectrum 
of 

^ERA ^ 1 (^MERA ^ pMERA ^ ^MERA ^ pMERA) ^ (^^gg) 



where descending maps V^f-^ and V^f-^ are those defined in (5.35). This 
tells us that homogeneous MERA states are also critical. 

It is commonly believed that, since MERA have larger stable causal cone 
sizes, they are more suitable to represent a state where the physics at very 
short ranges is sensibly different from the critical behavior at mid-to-long 



ranges: indeed homogeneous TTN states are forced to hold (5.66) at every 
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Figure 5.2: Critical exponents for an infinite XXY chain, calculated with 
a homogeneous MERA algorithm, as a function of the anisotropy A. Here a 
bondlink D = A was used. The red dashed lines report the theoretical values. 
Errors on the ground state energy (per site) appear in the inset. The graphic 
was kindly contributed by M. Rizzi and S. Montangero |29j . 



lenghtscale, even when £ ~ 2, while for a MERA state density matrices up 
to size 4 are loosely related related to their long-range physics. 



5.8 The importance of translational 
fluctuations in TTN 

Addressing translational problems with a non-translational variational ansatz 
might sound a sub-optimal choice for both analytical study and numerical 
simulations. Yet, TTN and MERA, as we just showed, manifest natural 
scaling properties that suit so smoothly critical systems, that are excellent 
variational candidates. Therefore, we are encouraged to wonder whether the 
capability of hierarchical Tensor Networks to reproduce a complete conformal 
symmetry, when L — )■ oo, goes beyond the mere paradigm of incoherent 
translational mixture. In other words: it is possible that in the limit of infinite 
layers fi (or sites), a tree Tensor Network recovers the same translational 
invariance that is forbidden to achieve at finite sizes? We dealt with this 
question by considering how the translational fluctuations renormalize and 
scale in tree networks, and developed a peculiar conclusion: fluctuations are 
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fundamental in TTN: they are necessary to describe interesting strongly- 
correlated physics, especially in the thermodynamical limit. If fluctuation 
vanish, the TTN state becomes trivial, and separable. We will now sketch a 
derivation of this relation between fluctuations and entanglement. 

In infinite systems, translational invariance is typically addressed as a 
hierarchy of equations, each one of them referring to a characteristic size i, 
and stating that the density matrices of that size are homogeneous in the 
lattice. Obviously, this implies that the same relation holds for any smaller 
size i' < i, thus the hierarchical relationship. 

Let us start from a one-site observable 0, and consider the translational 
fluctuation upon a finite ID lattice 



(5.69) 



where L = 2^^ is the system size. When we probe such fluctuation upon a 
homogeneous binary TTN state, the expression becomes 



Tr 



Tr 



e^e- (0" ( r/{,°^ ) -v^®v^ ( pf 



(5.70) 



where we used the fact that Tr[X]^ = Tr[X (g) X], while rj'^'^ is the same 



defined in (5.56). Moving to the thermodynamical limit is trivial now, and 
it reads 



ABl = lim (Ael'^l)' = Tr [0 ® e • (jT- pf ® pf )] . 



(5.71) 



with and pf being respectively the fixed point of ^ and V, which are 
uniquely defined once we assume that both maps ^ and V are mixing. We 
can now state that the thermodynamical TTN state is 'size-1 translational' 
if the quantity we just calculated vanishes; and it can be shown that this 
happens only if the density matrices and pf ® pf coincide. 
Indeed, let A and B be any two one-site observables, then 

Tr [(A ® S + S ® A) ■ (^r- pf ® pf )] = 

= AiA + B)l-AAl-ABl = 0, (5.72) 



which must be zero as as we are requiring that every translational fluctuation 
is negligible. Now let F be the swap operator F\a) ® = \(3) ® 
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F = F-^ = Ft. Clearly F{A ® B)F = B (g) A and Fipf" ® p^)F = ^ 
p^. Moreover is left invariant under the action of the swap gate, i.e. 
Fp^{FAF)F = p^{A). This implies that FjtIF is a fixed point of 
but since it is a mixing map the fixed point must be unique, which leads to 
j^l = Fjj'lF. The previous manipulations with the swap operator allow us 
to write the following equivalence 

Ti[A®B-{j^- p'^ ® pf)] = Ti[B®A-{j^-p'^® )] . (5.73) 



and since the sum of these terms is zero by (5.72), they must be both zero 
separately, for any operator A and B. But since tensor product operators 
generate the whole algebra of 2-site operators we must conclude that 

^=pr®Pr, (5.74) 

so that '0y and ® Vy must have the same fixed point. 

When this condition is verified then is automatically the fixed point 
of of and P/j as well, because 



= {p-V®V) (p?° ®p^^ 

= \{Vl- Vn){pT) ® {Vl - Vn)ipT). (5.75) 

Therefore, a sufficient and necessary condition for translational invariance to 
hold at size 1, is that Vl{pT) = Vr{p^) = V^p"^) = p^ . 

So far, so good. Let us proceed further and require translational invari- 
ance at size 2. The same derivation can be applied to two-adjacent sites 
observables Bj^+i, its result is that a common fixed point must be shared by 
maps 'D2,R and 1*2, l, where 

"J . , (5.76) 



It is clear that such fixed point must coincide with p^, since by (5.49) we 
know that I?2^2 = \{T^2,r + 1^2,l), and if we require that P2-s>2 is mixing its 
unique fixed point pf must be the shared fixed point of P2,_r and 1^2, l- At 



the same time, due to (5.75), the fixed point of T>2,r is pf (g) pf, and by the 



common fixed point property we must have that 

s{pT) = pT®pT = i?2^2(pr ® pT) = pT. (5.77) 



meaning that the two-sites reduced density matrix is a separable state, ac- 
tually a tensor product state. By extension, every reduced density matrix is 
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separable, and the system can manifest no quantum correlations at all, since 
Act defined in equation (5.62) would be the null operator. 

Resuming this whole discussion, we proved that by just requiring that 
the thermodynamical TTN-state manifests translational invariance property 
at size 2, we automatically end with a completely factorized state, all entan- 
glement is broken down. In this framework, we could state that translational 
fluctuations are needed in a tree network if we want to describe strongly 
correlated physics; they are unavoidable, even in the thermodynamical limit. 

It is curious to realize that the previous demonstration does not hold 
for a MERA geometry: this could actually be one of the first remarkable 
arguments for preferring MERA to Trees. Indeed for a MERA topology, 
according to this theoretical picture, we could hope to reproduce accurately 
conformal symmetry without implications of triviality for thermodynamical 
limit entanglement. 



5.9 Parent Hamiltonians of TTN states 

Finding ground states of physically-meaningful Hamiltonians has been for 
decades a maximal interest topic, it is also the very purpose of this whole 
Tensor Network variational ansatz itself. But even dealing with the reverse 
problem can be challenging and fruitful: given a quantum (many-body) state 
I^P), can we identify, characterize of even build a non-trivial Hamiltonian for 
which is the ground state? Of course, the research of such a parent 
Hamiltonian H, must be addressed in accordance to some physically sensible 
constraint: we might for instance require for H to be short-ranged, to be 
translational, or maybe capable of coupling only a limited number of particles 
per single interaction term. The more are the requirements, the harder is the 
problem. In my research work, I focused on Tree Tensor Network as ID many 
body states, and analyzed how to explicitly build a non-trivial, short-ranged 
and translational Hamiltonian which is parent for the TTN state, in PBC. 
In this section we will sketch the construction. 

Let us start from the formal definition of the Hamiltonian H we want to 
achieve, whose elementary terms H have limited size support i> <^ L: 

L 

H = Y, ^^0+1,... (5.78) 

In order for this Hamiltonian to be parent for our homogeneous-TTN state 
|vi>M(A)), we must ensure that its expectation values coincides with the min- 
imum of the spectrum of 7{, which is also the variational minimum of the 
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expectation values {H) on the whole Hilbert space of states: 

(*M(A)|?^|*M(A)) = min {^^^} • (5-79) 

At the same time, we could read H as an unnormalized translational aver- 
age. This allows us to summon again the incoherent translational mixture 
formahsm for TTN density matrices, which reads 

L 

|^-pll+Wo+J=2''IV[i/.p-M]. (5.80) 



=1 

L 



€o=l 

Now, suppose that exists a finite small (= non-scaling) u for which has 
non-maximal rank. If that is true, then has some nontrivial kernel, with 
strictly positive dimension, completely generated by an orthogonal set of 
vectors \kw)- Then we say that the elementary Hamiltonian term H is built 
as follows 

H ^^uj^\k^){k^\, (5.81) 

w 

with arbitrary positive weights cUyj > 0. With this prescription, wc obtain a 
Hamiltonian "H which is positive, as it is the sum of positive terms, nontrivial, 
if at least one cuu, is strictly greater than zero, and for which it holds 

(*M(A)|7{|*M(A)) = 2'^J2^^ MpI'^K) = 0, (5.82) 

w 

since every \Kyj) is in the kernel of pl^*'. But since is necessarily the minimum 
of the spectrum of as it must be a positive operator, |^['^](A)) is clearly 
a ground state for H. This is the idea of our construction. 

The central point of our proof, therefore, now becomes to demonstrate 
that for the (homogeneous) TTN state there is always some finite non-scaling 
size u for which the averaged z/-sites density matrix plf^ has non-full rank. To 
show this we will exploit the fact that descending maps can grow the size of 
density matrices, while the corresponding increase in entanglement is well- 
kept under control by isometricity condition on A (and X as well, in MERA). 
We will now discuss dimensionality relations leading to characterization of 
the smallest size u for which pl^^ is necessarily non-full rank. This result will 
obviously depend on the Tree or MERA geometry we are employing. In the 
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following instances, we will consider Trees and MERA having a renormaliza- 
tion dimension D which is equal to the physical local dimension d so that 
TN-homogeneity is meaningful up to the physical lattice. 

Binary Tree - in tree geometries it is not possible to establish any 
bound upon the entanglement of p\^\ for 2 is a (meta-) stable causal cone 



width. Let us move to size three: by adopting the formalism (5.42) we know 
that 



pr 



2^3 



1 



X®S 



V 



R 



P2 



(5.83) 



where X is the identical map, taking every operator into itself. Notice that 
even though the mapping S increases size, as it is the application of an isom- 
etry, it is left-invertible, and thus preserves the rank; by natural extension 
S ®X preserves the rank just as well. Now, the expression within the first 
square parentheses of (5.83) is a c?^ x (P matrix, whose rank is obviously 



bound by the amount of rows or columns (whichever the smallest), i.e. (P. 
Finally, it is clear that the maximal rank of a sum of two matrices A + B can- 
not overcome the sum of ranks of A and B separately. These considerations 
tell us that 



Rnk 



pT 



<2d' 



(5.84) 



where ^3'^' is a x matrix, and thus its rank is non-maximal whenever 
> 2d'^, which happens for a local physical dimension of d > 3. 
When we are considering a ring of 2-level systems, e.g. a spin-| chain, 

size three is not enough for ensuring non-maximality of the density matrix 

rank. Let us move to size four then: 



Pa -^^R 



5 



x®s®x 



1 



1 



s®s 



VR®Xi 



T>l[pV^ 



S®S 



P2 



(5.85) 



As previously motivated, both X ® S ® X and S ® S maps preserve the 
rank; and of course the ranks of their respective arguments is bound by their 
row/column dimension. Which leads to 



Rnk 



PI 



<d^ + d^ 



(5.86) 



which is always strictly smaller that d'^ for any nontrivial dimension d > 2. 
In conclusion, for binary trees, it is always possible to write a nontrivial. 
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translational, short-range Hamiltonian, according to the prescription (5.81) 



where every term H involves at most four (adjacent) sites, which become 
three for a local dimension d greater than 2. 

Higher branching TTN - when the Tree geometry has a branching 
number h higher than 2, the minimal size y of p"u \ for which it is possible to 
ensure non-maximality of the rank, grows. As an example, we can mention 
that for a ternary tree this size z/ is 5. But the growth rate, as a function of 
6, is somehow irregular, erratic; sometime the local dimension d is influent, 
sometime it is not. Nevertheless, for each geometry, a finite v definitely 
exists, and it is always equal or less than 26. In fact for a 6-branching tree 
network it holds 



Rnk 



P2b 



< h d'^' 



(5.87) 



which is always less than its row/column dimension ci^*. This interaction- 
range bound is not optimal, but it still does not scale with system size. 

MERA - The presence of the disentanglers also increases the minimal 
size of non-full rank pl^' . Precisely, for a binary MERA we get 



Rnk 



<2d* 



and 



Rnk 



< + d^- 



{5.i 



which tells us that we have to accept an interaction range = 5 for a local 
dimension d > 3, and move to u = 6 otherwise. If our MERA geometry is 
ternary, we have to push to seven sites, since 



Rnk 



pf 



< Sd". 



(5.89) 



At any rate, the effective range of our nontrivial translational parent Hamil- 



tonian might not be so short, but in the end the construction (5.81 ) is always 
possible in practice. 



5.9.1 Unfrustration and degeneracy 



It is meaningful to point out some general properties of the parent Hamilto- 
nians generated with the protocol (5.81) just described. First, we would like 



to highlight that these Hamiltonians are necessarily frustration-free, meaning 
that the TTN is ground state of every single interaction term. 



Indeed, let us consider the expansion (5.78). We built every single H term 



to be a positive operator so that the full Hamiltonian "H would be positive 
as well. But since the TTN state has zero expectation value of "H, it must be 



= (^M(A)|'H|^[^](A)) 



L 

4=1 



(vI/['^](A)|iJ,, 



+1, 



|*[^](A)). (5.90) 
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DOS 




Figure 5.3: Unnormalized density of states (DOS) of a TTN parent Hamil- 
tonian Us for on a spin-^ PBC chain of 8 sites [27j. Here we used the sample 
isometryA : |0) |01), |1) 2~ ^/^(|00) + |ll)). The interaction term if was 
calculated from pf{A), via (5.81), with random positive weights u^. Notice 



that the ground space is 32-fold degenerate: its dimension is twice the lower 



bound given by (5.96), the x2 factor deriving by a hidden Z2 symmetry. 



Now, the only way for a sum of positive terms to be zero is that every term 
is separately zero: 

(*M(A)|ii,„+i,...,,„+,|xl>['^](A)) = V£o, (5.91) 



i.e. "H is unfrustrated. This one consideration is curiously related to refs. [621 
[63] , where it was shown that for structured frustration-free Hamiltonians, it 
is possible to build analytically a ground state via Tensor Network designs. 
In some sense, our result is the other face of the same coin. 

Another manifest property of our TTN parent Hamiltonian, and some- 
how related to the unfrustration, is that the ground space of is highly 
degenerate, and we can characterize it to some extent. To show this, we will 
consider for simplicity the case of a Binary Tree and a local dimension li > 3, 



and assume that the construction (5.81) is being developed at the thermo 



dynamic limit /i = 00. We will also request, as additional hypothesis, that 
has full rank. This assumption is typically weak, practically guaranteed 
in numerical settings, where stochastic noise makes every matrix full rank. 



Now, via (5.81) we derive a positive interaction term H having support 
in the kernel of pf. To begin with, let us prove that H is in the Kernel of the 
ascending map ^3->2, adjoint to '^2^3 of eq. ( |5.83[ ): by definition we have 



= Tt[H- p^] = Tt[H- ©2^3 (p^)] = Tr [^3^2 (H) ■ p^] . (5.92) 
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where since ^3_j>2 map is completely positive and unital, ^3^2 (H) is surely 
a positive operator, is also positive, and being full-rank, its smallest 
eigenvalue Aq is strictly greater than zero; then 

= Tr [^3^2 (H) ■ p^] > Ao Tr [^3^2 (H)] > 0. (5.93) 

So it must be that the trace of ^3_s.2 (H) is zero, but the only positive trace- 
less operator is the null operator, thus ^3-s.2 (H) = 0. This very argument 
will let us characterize ground spaces of Hamiltonians H generated by this 
interaction term H. 

In fact, let us consider a finite system now, with the same local dimension 
d and an ever number of sites L = 2i. we will define our trial state as a generic 
pure state on i sites l^^o), that we grow to 2i sites by means of a single layer 
of isometries A: the same A tensor we used to build H. Let us sketch the 
Tensor Network design of this trial state as follows: 




(5.94) 

It is trivial to show that this trial state is definitely a ground state for the 
Hamiltonian Ti = J2e =1 -^^o+i---^o+3- Indeed, let us write 



{n) = Tt[h 



pi] 
Pi) 



pI] 



0, 



(5.95) 



where pi (resp. pj,) is the reduced density matrix, i/-sites translationally 
averaged, of the trial state before (after) applying the layer of isometries A. 



Equation (5.95) is telling that the trial state is a ground state of the Hamil- 



tonian, regardless from |\E'o)- Actually every state that is written in the form 



(5.94) is a ground state of the system, and the layer of isometries preserves 



orthonormality, so we can identify (at least) a set of d orthogonal ground 
states. In conclusion, the ground space of any "H has a wide degeneracy, 
namely: 

dim [Ker (H)] > (5.96) 

with L the size of the system. The same discussion can be applied to the 
MERA case, and leads to the same result, although the ground states this 
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time are built by attaching to |\l/o) a full MERA layer, with renormalizators 
A and disentanglers X together. 



5.10 Open boundary MERA 



In section |572l we mentioned that Tree Tensor Networks are equally suitable 
to simulate open boundary systems as well as periodic boundary systems. 
This is mainly due to the fact that there are pairs of adjacent sites that 
renormalize separately for an arbitrary number of layers, thus actually con- 
stituting an inner-boundary when approaching TD-limit (honestly, a thermo- 
dynamical TTN is geometrically equivalent to the frontier of the Cantor set). 
At the same time, we stated that MERA have a natural attitude for periodic 
topologies, since a single MERA layer couples every pair of neighboring sites. 

It is easy to see, however, that a MERA state needs only a little adjust- 
ment to its network geometry, to take properly into account the presence 
of an open boundary. This is naturally done by embedding the boundaries, 
interpreted in the most general setting as a pair of additional degrees of 
freedom (ancillae), in the MERA picture, while preserving the bulk network 
pattern, and thus the critical bulk properties. The present idea leads to the 
following network design: 




(5.97) 

which we refer to as (open) boundary-MERA [28]. As you see from the 
diagram, during the same algorithm-step when we apply disentanglers, we 
also allow the boundary ancilla to couple locally with the system. We will 
add the constraint that this coupling operation is represented by a unitary 
gate. Such additional requirement is necessary to preserve the causal cone 



relations we argumented in section 5.3 Clearly, isometricity and unitarity 
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conditions still hold for renormalizers A and disentanglers X as before, i.e. 





(5.98) 

where the bottom tensors are complex conjugate and up-down reversed ver- 
sions of the ones standing directly above. Preservation of the causal cone 
property, together with the fact that the binary MERA pattern is identical 



to the PBC case in the bulk, tell us that descending map equations (5.35) 



still hold when the three involved sites (causal cone width of 3 is still stable) 
are far from the boundaries, i.e. 



[tA 7-) f 

P2^+l, 2^+2, 2^+3 ~ ^3'^3,R yPe/+l/+2 J 



and 



(5.99) 



where maps T>^^^ are those of eq. ( 5.36 ) . When approaching the boundaries 



we must accordingly define descending quantum channels that involve the 
new system-ancilla coupling element. For instance, close to the left boundary 
we get: 



Pl,2,3 



/c, 



Paa,2 



and 



Pa,1,2 



Br 



Pa,1,2 



(5.100) 



where the subscript A refers to the left-ancilla degree of freedom, and the 
completely positive trace-preserving maps JCl and Bl are given by 



L 




B 



L 




(5.101) 

Similarly, the maps ICr and Br, ruling the layer-recursive relations of density 



matrices at the right boundary, are the left-right specular versions of (5.101 ). 
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Notice from (5.100) that once a causal cone touches a boundary, it sticks to 



it along its upward propagation. 



5.10.1 Consistency of the thermo dynamical limit 

Before speculating on how the presence of the boundary influences the sys- 
tem, we want to show that as we approach the thermodynamical limit (and 
consider translationally-averaged quantities), we recover the same physics 
of the corresponding PBC-MERA. In particular, we will assume that the 
boundary-MERA is homogeneous, i.e. renormalizers A are identical, and so 
disentanglers X, and even boundary- ancilla couplings; then move to /i — )■ oo. 

We will sketch the consistency proof for size 3, which extends trivially to 
any size thanks to the shrinking properties of the causal cone width. Consider 
the translationally averaged 3-site reduced density matrix 



1 



4=1 



^4,4+1,4+2- 



(5.102) 



We are actually not considering the ancillae as part of our system in this 
framework, but the result would not change even if we did. We can now 



adopt recursive relations (5.99) and (5.100) and write: 



Pa,1,2 



Pl-i,l,A' 
1 



+ ^1 - ^^irfn j ■ ^3.3 (p-r' j . (5.103) 

In order to identify the thermodynamical limit of this quantity, it is important 
to notice that the following expression converges to zero in trace-norm as fi 
grows: 

4 



lim 



Ps 



Ps 



< lim 



2M+2 



0, 



(5.104) 



pT 



lim pg 



where we used the fact that CPT maps are contractive, that the trace norm 
of a density matrix is 1, and triangular inequality. The previous equation 
guarantees that pg^' goes in the limit to the fixed point of T's-j.s, as 

^[^1 - lim V,^s (p^'O = ^^3^3 (pT) , (5.105) 

but the fixed point of 'Ps-j.s also characterized the thermodynamical limit in 
the periodic MERA, and since we assume for 2^3^.3 to be mixing, it has a 
unique fixed point, thus the two states must be the same. This proves that 
PEG and OBC MERA manifest the same averaged physics in the thermody- 
namical limit, which is the TD-consistency argument we requested. 
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5.10.2 Boundary fluctuations and permeation 



Boundary conformal field theory prescribes a direct relation between criti- 
cal exponents of two-point correlation functions in the bulk, and one-point 
fluctuations close to the boundary of an OBC critical system. We are now 
going to investigate such one-point expectation values as a function of the 
distance from a boundary, say the left one, and compare these analytical 
results with our previous acknowledgements involving correlation functions 



in MERA, namely (5.66) and (5.67). 



To do this, we consider the expectation value of a 3-site observable ap- 
plying at a distance £ ^ L ~ oo from the, say left, boundary: 



lim (6 

fl—^OO 



£,i+l,i+2/ 



lim Tr 



/I— >oo 



e- 



(5.106) 



which appears as a function, upon £, on how the influence of the boundary 
permeates inside the (infinite) system. We will now assume, for simplicity, 
that such distance £ is a power of 2, say £ = 2^^. Then, by adopting the 



formalism of (5.99), we obtain 



<P~2,(e) = lim Tr 



^ ^3^3,L lPl,2,3 



lim Tr 



(5.107) 



where p'^^ 2 is the fixed point of Bl, unique if the map is mixing. The only 
residual dependence on r is in the number of times the map A^^s^l is to 
be applied to the operator (or, equivalently, V^^^ l is to be applied to the 
boundary density matrix). We can now proceed by adopting an argument 



similar to (5.65), i.e. let us assume that Qa,L is an eigenoperator of ^3^3,^, 



then the permeation function obeys an exact power-law behavior: 



(5.10J 



where is the relative eigenvalue, |Aa| < 1. It is evident that permeation 
functions near the boundary are critical indeed, as the ancillary degree of 
freedom strongly correlates with the system. By expanding a generic observ- 
able 9 in the generalized eigenoperator basis of ^3^3,^ we can write 



[As- 



'(l0g2^), 



(5.109) 



d.w 



with critical exponents C« determined by the spectrum of ^3-5.3,^ via loga- 
rithmic relation C,a = log2 Aq. 
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Now we would require somehow that the previous result does not hold 
only for 2^ distances, but that it generalizes smoothly to other distance, 
so that our boundary-MERA ansatz would be able to capture conformal 
symmetry as closest as possible. In particular, let us now consider distances 



of the form i = 2"^ — 1. Then, (5.99) tells us that the correct recursive relation 
reads: 

P2^-l,2^,2^+l = ^3^3,R (p2^-iLi,2^-i,2^-i+l) ' (5.110) 

Notice that now is the other size-3 descending map (Vs^s^r and not Ps^s.l) 
that rules the recursion, unlike ( |5.107[ ). Thus, the permeation function be- 
comes 

W=2^-iiQ) = Tr [AlZ^l^ (6) ■ ICl (pT,i,2)] , (5.111) 

stating that now the critical exponents are logarithms of the spectrum of 
A3^3,R, and its eigenoperators Qa,R are one-point primary fields. 

Let us introduce the translational regularity requirement. For instance 
we demand that critical exponents and primary fields do not depend on the 
point i at which the permeation function is computed. It is easy to see that, 
when accepting this, it immediately follows that 

P3^3,L = P3^3,R = ^^3-^3, thuS 1 1 2l 



But now, recall that by (5.66), the ^ was the map ruling the two-point 



correlation functions in the bulk, two point critical exponents given by 



the logarithms of its spectrum. Thanks to (5.112) we know that two point 
primary fields in the bulk correspond to the one-point ones at the boundary, 
and critical exponents satisfy 

Ca = \ ^a. (5.113) 

As was pointed out by P. Calabrese [351 EH], this is one of the fundamental 
properties prescribed by boundary conformal field theory. One-point critical 
exponents at the boundary are exactly half of the two-point ones in the 
bulk. We re-derived independently this feature by solely exploiting geomet- 
rical features of boundary-MERA network and smoothness requirements. 



5.11 Hybrid MPS ^ TTN networks 

When a quantum many-body system draws near to a second-order quantum 
phase transition, noncritical and critical properties begin to overlap. Often, 



144 



CHAPTER 5. TREES AND MERA 



as system parameters approach the critical region (especially in ID settings, 
where phase transitions are allowed only at zero temperature) islands char- 
acterized by strong correlation within start to appear, out of a noncritical 
long-range landscape. Then the characteristic size of these regions them- 
selves increases, until they reach the lenghtscale of the whole system when 
the critical point is achieved. 

In order to simulate efficiently this type of quantum behavior with a 
Tensor Network ansatz, one would like reproduce both power-law like corre- 
lation scalings up to a tunable finite distance, and exponential decay rates 
beyond. A promising candidate for this goal would be a structure which em- 
beds the self-similar geometry of Matrix Product States, which express non- 
critical character in a natural way as we saw in section 3^, when observed 
at large distances, and resembles a hierarchical Tensor Network, either TTN 
or MERA which both bear strong correlation capabilities, when observed in 
proximity of the physical bondlinks. Such idea leads, almost obviously, to 
the design of a Hybrid MPS^TTN structure, picted as follows: 




(5.114) 



where the system size L is given by the interplay of two parameters: n, the 
number of Matrix Product blocks, and fi, the full depth of the binary tree 
^ curtain^ hanging from the MPS layer. So that L = 2'^n. In order to suit as 
best as possible translational invariance, despite the presence of Tl-breaking 
tree geometry, we will consider the MPS tensors A to be homogeneous, and 
the isometrics A of the tree as well. It is clear that, as we want for the system 
to reach thermodynamical limit but keep strongly-correlated islands at finite 
size, we will fix the curtain height /i, and increase n — )■ oo, so that L goes to 
infinity too. 



The algebraic ingredients in our analytical study of (5.114) hybrid net- 
works are again the translationally-averaged density matrices p^i \ a[^'^ and 



r]f , respectively defined in (5.37) and (5.56). Recursive relations involving 



these elements will let us figure out correlation behavior and scaling-laws. 
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MPS as starting point - In section 3.4 we learned that the fastest 
way to keep under control normalization and physical quantities in thermo- 
dynamical limit MPS is to choose the appropriate peripheral gauge, say the 
left one. When fi = 0, the involved density matrices are completely defined 
by the matrix product properties. Namely 



pfi=^(<|.+|^®A;^|Q)|.)( 

r,s=l 
d 



rnl ^ (5.115) 

crf= Yl i^^\{As,®A:.^)E'-'{A,,(g^AX)\Q)\s,S2){nr,\ 



ri,r2,si,S2 



where the MPS matrices {^4^} are chosen so that they satisfy the gauge 
condition A-l^s = 1. and E = ® A* is the transfer matrix of the 

Identity operator. Correlation vectors \Q) and ($+| are respectively the right 
and left fixed points of E, they are both positive when read as matrices, and 
in particular ($+| = J2a(^^\ corresponds to the identity itself. 

Adding the TTN curtain - we can now increase /i to a finite nonzero 
value. Recursive relations between density matrices standing different tree 
layers are as usual given by descending CPT-map formalism. So that 

pf =i)(pr]) 

aM=^(af-^l) (5.116) 
for every < u < /j,. It is clear that, while p^i^ can be always expressed in 



terms of (5.115) as V'^{p^"^), expressions for cr|^^ and rj^^ depend whether the 



distance i is smaller or larger than the island size 2^. 

Small scale regime {£ < 2^) - For simplicity, let us choose i = 2'^, 
clearly with p < fi {p integer). Then the appropriate recursive relation read 

We still have to exhibit an expression of a^i^ and ri^^\ for < u < fi. To do 
this in a clever way we introduce a new density matrix 9^'^\ defined as 



2Pn 



4=1 
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Interestingly enough, both cr[ and r/^ can be extracted from 9^"^ by partial 
trace, precisely 

= Tr3,4 [6^"^ and t]^^^ = Tr2,3 [6^"^ ■ (5.119) 
At the same time, ^'^^ satisfies the following recursion 
,M = w ^ + , (5,120) 

starting at the MPS layer from 6^^^ = af^ (g) af'K 

Large scale regime {i > 2^) - Things are easier now, provided that we 
choose, again for simplicity, i = (r integer), with no finite-size effects 
£ <^ L ~ 00. The correct expressions are just as follows: 

aS=^n4°'), (5.121) 
vi:l=0'iv^^^)- (5.122) 

We formally achieved every quantity we are interested in. We are now ready 
to combine together the previous results to investigate correlation properties 
of Hybrid Tensor Networks. 



5.11.1 Two-point correlation functions 

We are interested in calculating two-point correlators at fixed distance and 
translationally averaged on variational states of the form (5.114). Let us then 
recall 



1 ^ 

l?ie, 0') ^ - ® ©k+^i) - (0N)(0k+^+i]) 



Tr 



(5.123) 



the appropriate correlation function of and 0' at distance i. Let us focus 
on the scaling laws of this quantity, in either of the two regimes we mentioned, 
rather than an exact formal expression for it. 

Small scale regime - Again i = 2^, with p < ^. By using the recursive 



relation (5.120) for 6^'^^ we can write 



4=2.(0>0')=Tr 



(0 ® 0') o (Tr3,4 - Tr2,3) [W^-^ (^M)]] . (5.124) 
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the residual dependence on p (and fi) is left only in the number of times the 
maps and W are to be applied. This allows us to exploit the Jordan block 
expansion scheme for multiple actions of the same matrix, as we did before, 
leading us to 



=2P 



^{A}(P) 



(5.125) 



Aw 



where (resp. Aw) are eigenvalues of ^ (of W), and V^x} are finite-degree 
polynomials. Now recalling that p is the logarithm of the distance i, it 
is clear that (5.125) is dominated by power-law behaving functions, apart 



logarithmic corrections, whose (quasi-critical) exponents are determined by 



Xo 



logs 



Aw 



(5.126) 



Large scale regime - The system behave differently when i = 2^r, for 
1 < r ^ L. In fact at these distances eq. (5.123) becomes 



^^^(e, e') = Tr [(6 ® 0') -0^ {al 



Vr 



E-|g). (5.127) 



In the end we are actually calculating a r-fixed distance correlator upon the 
MPS layer of the observable X''(0 ® ©')• But since we chose a homogeneous, 
left-gauge MPS description in the thermodynamical limit, we can perfectly 
recover the corresponding result ( 3.28[ ), 



I.e. 



£2..r(e,e') 



|Ae|<i 



(5.128) 



In fact, instead of (5.125), now r is proportional to i and we obtain an 



explicit exponential decay, going to zero at infinite distance since the sum 
spans only the eigenvalues Ae of smaller of 1 in modulus. 

In conclusion, we discovered that the hybrid Tensor Network state geome- 



try, designed in (5.114), manifests a quasi-critical character. At short ranges. 



strong correlations identified by power-law two-point functions, arise; their 
exponents characterized by the isometry element A of the TTN portion. At 
long ranges, the behavior is evidently noncritical, correlations vanish expo- 
nentially and ruled by the MPS block A. The number of TTN layers we 
adopted in the scheme effectively determines the lenghtscale io of the strong- 
correlation islands £o ~ 2^. 

MPS ^ MERA Hybrids - The previous discussions were formulated 
for a binary tree curtain attached to the MPS-basis layer. But it is easy to 
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see that the same results can be generahzed to other tree geometries and for 
MERA as well. So even if we had a network structure of the form 




we would similarly obtain a quasi critical regime. In fact, descending maps 
for MERA replace the role of those for TTN, and the scaling laws do not 
sensibly change. 



5.12 Higher dimensions and comparison 
between TTN/MERA and PEPS 

So far in this chapter, we discussed on how tree Tensor Networks and MERA 
can be used as a suitable and efficient variational for many-body ground 
states of one dimensional system. But it is clear that their power is not 
limited to ID systems only, indeed generalizations of Trees and MERA to 
higher dimensions is natural and intuitive. 

The idea remains the same [61]: we may think to perform a real-space 
renormalization group of the quantum lattice state, expressed in the language 
of density matrices, for any amount of spatial dimensions t^D. By perform- 
ing this iterated process, we obtain a class of states which is equivalently 
expressed as a Tensor Network, namely, a tree network in #D dimensions. 
For instance, let us consider a 2D square lattice, then the simplest 2D tree 
geometry we can think of is the one that maps a 4- adjacent sites plaquette 
into a single renormalized site (a). The renormalization element A becomes 
then a 5-link tensor, and read as a D'^ x D matrix is an isometry A"''A = 1. 

Alternatively, again dealing with a square lattice, one could prefer to 
renormalize together pairs of sites which are horizontally adjacent (so that 
we are actually coarse-graining only one of the two dimensions), and at the 
following step those who are vertically adjacent (b). Clearly, the latter tree 
geometry has the advantage on the former that is easier to contract and 
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thus more efficient, for the same refinement parameter D. At the same 
time description capabihties of the second choice (in terms of variational 
manifold) are reduced. In diagrammatic expression these two TTN designs 
read respectively: 





(5.130) 



The possibilities are many. The user is encouraged to adopt for her sim- 
ulation the T^D TREE/MERA geometry that suits most symmetries (see 
appendix |B| and other global properties of the problem under study; for 
example, a plaquette-tree geometry would embed more naturally the square 
group symmetry. 

Many of the algebraic properties we discussed in this chapter for ID 
TTN/MERA are preserved by their higher-dimensionality versions. The 
causal cone property is one of these: as TTN remain Tensor Networks with- 
out closed loops, transforming them according to the peripheral gauge (with 
respect to the hat as nucleus) is always possible without loss of generality, no 
matter the geometry or dimensionality. Thanks to the causal cones, it is clear 
that even in high-D, the contraction of a Tree network, or MERA, stays effi- 
cient: the number of required elementary operations scaling logarithmically 
with the volume. 

This is probably the most convincing argument for preferring TTN/MERA 
to PEPS, from 2D up. Despite the two families of Tensor Network mani- 
fest promising description capabilities for physically-relevant ground states, 
PEPS are computationally complex objects, while hierarchical Tensor Net- 
works keep their high efficiency. It is worth mentioning that it was shown 
in [17| that starting from 2D up, MERA (and thus Trees) satisfy an en- 
tanglement area law, and therefore can be efficiently mapped, through a 
well-defined formal algorithm, into PEPS. 
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In the end, Trees and MERA are simulation tools hard to ignore, for 
any physical setting: in ID they are the answer on how to simulate critical 
systems in a natural and cheap way; in two or more spatial dimensions, 
they constitute an intriguing sub-class of finitely correlated states that are 
efficiently contractible. 



Conclusions 



In this thesis we introduced, developed, and analyzed a large class of vari- 
ational tailored quantum many-body wavefunctions, called Tensor Network 
states. These states are meant to be used as variational ansatze for interact- 
ing particles problems on a lattice; minimization algorithms adopting these 
trial functions do not require any a priori knowledge on the model we are 
studying. Predicting the success of Tensor Network states for a certain set- 
ting is based upon arguments borrowed from quantum information theory: 
primarily entanglement. 

We discussed profoundly how Tensor Network designs are somehow the 
variational counterpart of some numerical renormalization group procedure. 
In fact, the manifold of quantum many-body states that can be constructed 
via a RG-algorithm, can be equivalently identified by a tailored analytic ex- 
pression, where the variational descriptors are tied together through simple 
linear algebraic relations. For these reasons Tensor Network states preserve 
all the faithfulness and simulation power of numerical renormalization groups, 
yet the variational picture presents several practical advantages. Two im- 
portant example are: a more immediate way to access physical information, 
and a more free numerical manipulability with a consequent computational 
speed-up. 

In particular, we showed in chapter |2] that Matrix Product States arise 
from White's DMRG. Since Matrix Product States correspond to finitely- 
correlated states in ID, they are allowed to manifest only the correct amount 
of entanglement of ID non-critical ground states, and at the same time to be 
extremely efficient for computation. This explains why DMRG (and MPS- 
based algorithms alike) are so successful in one-dimensional systems. 

In chapter [3] we presented a generalization of MPS formalism that extends 
to Periodic boundary systems, with the care of keeping the right amount of 
allowed correlation. This led us to the definition of a homogeneous PBC-MPS 
representation, which was instrumental in the definition of thermodynamical 
limit for Matrix Product States. Once we defined an infinite MPS, we could 
investigate the correlations behavior of these states, and showed that under 
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the assumption of keeping a finite refinement parameter, the target state is 
distinctly non-critical. 

Matrix Product States are considered the fundamental template that led 
to the more general definition of Tensor Network states. Various classes of 
Tensor Networks manifest properties depending strictly on their geometry. 
Still, it is possible to determine general features common to every TN-state; 
those we sketched in chapter |4} 

Finally, chapter |5] contained most of the analytical advancements I con- 
tributed. There we discussed about Tree Tensor Networks, and MERA, the 
latter geometry presented as an enhancement of the former, but still sharing 
most of its scaling properties. We explained how these TN-variational states 
arise from a real-space renormalization group technique, an origin which 
leads unavoidably to scale invariant physics for these states. In fact, we 
showed that by adopting a CPT map formalism, we can well-define the ther- 
modynamical limit for TTN/MERA states, and this state manifest critical 
behavior, identified by power-law decay rates of two-point correlations. 
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Appendix A 



Completely Positive Trace 
preserving maps 



Here we list some features of Completely Positive Trace preserving (CPT) 
maps, as they are common and useful tools in quantum information theory, 
and thus thoroughly studied [65] . They represent quantum channels, having 
the property that they map density matrices into density matrices, and of- 
ten describe the time evolution of quantum states when both coherent and 
incoherent sources couple with the system. In Tensor Network settings, CPT 
maps are primarily used as inverse transformation of some numerical renor- 
malization procedure, typically applied to a density matrix, and thus work as 
propagators for renormalized density matrices within the network structure. 

A.l Definition 

A homomorphism between matrix spaces Ai : C"^" — ?■ C™^"* is said to be a 
CPT map if satisfies the following requirements 

1. M is positive: WA>0,Ae C"^" ^ M{A)>0 

2. Ai is completely positive, meaning that its identity extension on any 
larger space is still positive. Let us define Ai : C''^'' £nxn _^ 
Qfexfe ^ Qmxm defined on separable matrices [Ik (g) M]{B (g) A) = B (g) 
M{A). Then the complete positivity reads: VC > 0, C G C^^'' (g) 
£nxn ^ ^ M]{C) > 0, for every extension k. 

3. Ai preserves the trace: Tr[A^(A)] = Tr[74]. 

Of course, a matrix A is positive if its spectrum lies in IR^, or also {ip\A\ilj) > 
W\ip) if we are in a Hilbert space. Choi's theorem [6l] states that a mapping 
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between matrices Ai is completely positive iff it exists a set of Kraus oper- 
ators {Vs}s for Ai, namely a set of matrices Vg G C™"^" which satisfy: 

5^\/;v; = lnxn and ^K^V; = A^(^) VA (A.l) 

s s 

For later reference, it is also convenient to define also the map adjoint of A4 
with respect to the trace scalar product between matrices {A\B) = Ty[A^B]; 
by definition M^'^^ : C'"^'" ^ C"^" so that 

{A\M'''^^{B)) = {M{A)\B) = TT[{M{A)yB] = ^Tr[K^V;5]. (A.2) 

s 

Then, by using cyclicity of the trace, you immediately see that JH'^'^^ must 
read: 

M'^'^KB) = J2ylBV, V5. (A.3) 

s 

Which tells us that Ai^'^^ is again completely positive, and unital, i.e. it maps 



the identity operator onto itself, thanks to (A.l), but not necessarily trace 
preserving. 



A.2 Spectral properties 

When CPT maps Ai are endomorphisms {n = m), they can be expressed as 
square matrices, and expanded in a basis of eigenoperators and generalized 
eigenoperators as usual. However, spectrum and eigenmatrices of a CPT 
map are bound to undergo certain properties: 

• Every eigenoperator 0\ with eigenvalue A 7^ 1, must have null trace. 
This is clear from the fact that Tr[OA] = Tr[A^(OA)] = Tr[AOA] = 
ATr[OA], and since A 7^ 1 its only solution is Tr[OA] = 0. Also, it 
follows that such a 0\ can not be positive, because the only positive 
traceless matrix is the null matrix. The traceless requirement trivially 
extends to generalized eigenoperators as well. 

• The spectrum of M. must he symmetric with respect to the real axis. 
Indeed, M{A^) = Y^s^sAW} = T^s^^sAV})^ = (A^(A))t. Then, if 
0\ is eigenoperator with eigenvalue A, we have A^(0|) = (A^(Oa))^ = 
(AOa)^ = y0{, so also A* is in the spectrum. In particular, a hermitian 
eigenoperator has necessarily a real eigenvalue. This also tells us that 
the spectrum of M. and A^^'^J coincide. 
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1 always belongs to the spectrum of Ai. By absurd if Ai had no eigen- 
value 1, then its entire basis of generahzed eigenoperators would be 
traceless, but that would necessarily be an incomplete set because the 
identity operator 1 can not be generated. Another way to see this is 
that 1 is always eigenoperator of A^^'^j, and 1 its eigenvalue; but spectra 
of A^^'^j and Ai are equivalent, so 1 is also eigenvalue for Ai. 

The generalized eigenspace of X = 1 coincides with the strict eigenspace 
In fact if we assume that there is a Jordan block of dimension 2 
or greater, we can define the generalized eigenoperator 0[ for which 
M'^^iO'^) = 0[ + 1. but the equation a = Tr[0[] = Tt[0[] + Tr[l] = 
a + n has no solution. 

The spectral radius is 1, or equivalently, CPT maps are contractive. 
To prove this, we are going to show that if Ai^'^^ had an eigenvalue A 
greater in modulus than 1, it could not be a positive map. Indeed let 
Ox be the related eigenoperator, and consider 0+ = Ox + 0\, which is 
hermitian, but traceless and thus not positive, meaning that it exists 
a (normalized) vector lip) for which {ip\Oj^\ip) = w < 0. Now, the 
operator (1 + eO+) is definitely positive for < e < ||0+||. Therefore 
[M'"^f{l + eOx) = 1 + eiX'^Ox + A^^Oj) = Tg should be positive for 
every q, but 

{tfjlTgltfj) = 1 + £:| A|^ {w cos{(j)q) + z sm{(j)q)) , (A. 4) 

where (p = arg(A). This expression either oscillates with exponentially 
increasing amplitude, with period 27r/0, or it is monotonically, and 
exponentially, decreasing if = 0, since w < 0. Either way, sooner or 



later we encounter some integer q for which (A. 4) is negative, telling 
us A^^'^j mapped a positive operator into a non-positive one, which is 
the desired absurd. 

This final remark on contractivity of CPT maps is particularly impor- 
tant, it tells us that within the space of matrices, there is a proper subspace 
which is attractive, i.e. that every operator is driven towards it by multiple 
applications of the map Ai, its distance from such set exponentially decreas- 
ing. Also we know that the attraction subspace is generated by generalized 
eigenoperators of eigenvalues A of modulus |A| = 1. 

Obviously, the case that gathers greatest interest from this point of view 
is when this subspace is one-dimensional, so that the map contracts the 
projective space into a point. We will see that this additional request is also 
related to spectral properties of A^. 
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A. 3 The Mixing requirement 

A CPT map Ai is said to be mixing when infinite applications of A4'^ = 
JH o . . . o J\A contract tlie wliole space of matrices (modulus the trace) into 
a unique point A, i.e. 

lim M'^iA) = A Tt[A] VA. (A.5) 

It is clear that when this condition holds A is obviously a fixed point of the 
map, since 

M{A) = o lim M\l/n) = lim M''+\l/n) = A. (A.6) 

q—^oo >oo 

which immediately tells us that A is positive and Tr[A] = 1, so A is a density 
matrix. The mixing condition has consequences on the spectral periphery of 
Ai , which are also related to well-known physical properties [66] : 

• Ergodicity - The eigenvalue 1 of Ai is simple, meaning that the re- 
lated eigenspace is one dimensional, i.e. Ai has a unique fixed point. 
When this statement fails, it is clear that it exists a whole manifold 
(at least one dimensional) of operators On, with Tr[OM] = 1, which 



are all fixed points of the map. Thus eq. (A.5) breaks down since 
A^°°(Oh) = Om and Om not unique. 

• Relaxation - Ai has no modulus |Aq| = 1 eigenvalue Xa, except for 
1 itself. This tells us that every component of a given observable O, 
expanded in the generalized eigenbasis, other than A, decays exponen- 
tially with a rate governed by the second greatest modulus eigenvalue 
IA02I < IA2I < 1. By contradiction, if an eigenvalue A = e**^ ex- 
isted, one can construct a whole set of operators that rotate infinitely 
around the fixed point A, like 

Al^(A + aO^) = A + ae*^^0^ (A.7) 

and the sequence in q has no limit. 

When the map Ad is mixing, expanding the action of several, but finite, 
applications of Ai in the generalized eigenbasis Oa,d,w (« - eigenvalue index, 
d - Jordan block index, w - position within the block) is particularly useful, 
and it reads 

Ai'^iA) = ATt[A] + J2Kj2'^t'~'"^{q) O.A^, (A.8) 

a=2 d,w 

with 'P|^'(g) is a polynomial function of degree x, with coefficients depending 
on A and Aq,; obviously all the sum terms vanish at g — )■ 00. 



Appendix B 

Symmetries in Tensor Networks 



Exploiting symmetries in numerical simulations of quantum problems is one 
of the best and most valuable techniques we can adopt to drastically improve 
the efficiency of computation, with no actual loss of accuracy. When the 
Hamiltonian of the many-body problem is invariant under a group of unitary 
transformations, then an immediate characterization of the ground states, 
under the action of the group itself, emerges. Even in the rare cases when 
a discrete symmetry is spontaneously broken, it is always possible, thanks 
to the wave-like formulation of quantum mechanics, to identify a symmetry- 
invariant ground state. 

In variational contexts, where the efficiency of every protocol is extremely 
sensitive to the number of effective parameters, the total amount of such 
descriptors is heavily reduced every time a symmetry constraint is embedded 
in the framework. Of course, the capability of upholding a symmetry in a 
given variational ansatz is not guaranteed a priori: the manifold of states 
belonging to the variational family must be able to recognize, distinguish, 
and capture the physics of the symmetry group as a whole. 

Tensor Networks, as we discussed thoroughly, are variational counterparts 
of numerical renormalization groups. Since RG-processes naturally manifest 
a notion of locality preservation, TN-states seem most suitable to reproduce 
symmetries which are global, but act locally as an uncorrelated product of 
on-site transformations. The strategies to implant symmetries within Matrix 
Product States had been known from the DMRG era ||5], [67] ; generalizations 
to other Tensor Network geometries were studied in depth by S. Singh et 
al. |6H1 |69]. It was shown that the overall effect of inserting a symmetry 
constraint in a TN-ansatz is that one can easily isolate residual variational 
degrees of freedom out of structure factors and selection rules. In practice, 
every tensor decomposes, into a variational and a structural part (or frag- 
ment), which are again tied together by network linking relations. 



157 



158 



APPENDIX B. SYMMETRIES IN TENSOR NETWORKS 



B.l Point wise symmetries and 
representations 

We will now define the class of symmetry transformation groups we will take 
into account through the following sections, and then review some basic, but 
useful, results of representation theory. Let Q be our compact symmetry 
group; it can be a finite group as well as a Lie group: typical candidates are 
Z„, 0(n), SO{n), U(n), SU(n) but also dihedral groups are common. We 
identify {Ug}g^g as its unitary representation on the one-site (i-dimensioned 
Hilbert space: Ug^^g^ = Ug^ ■ Ug^ Then the overall transformations are given 
by tensor products of on-site unitaries 

L 

Uf'' = ®Uf\ (B.l) 
1=1 

where L is the total number of sites, and Ug does not depend on the site 
I on which it acts, as we are requiring homogeneity of the local representa- 
tion. In literature, these transformation are often referred to either as local 
symmetries, as U®^ is a product of local terms, or as global symmetries, as 
U®^ has global support. To avoid misunderstandings, we prefer to call trans- 



formations of the form (B.l), pointwise symmetries, as they are an overall 
operation acting on every site singularly. Pointwise symmetries are extremely 
relevant from a physical point of view: they can uphold most of the extensive 
constraints of a physical framework, like total particle conservation, parity, 
or spin conservation. Trivially, the f/®^ form a group and are an actual 
representation of Q, although it is hardly an irreducible representation, even 
when the original Ug is irreducible. 

Before moving further, let us briefly recall some basic principles of ma- 
tricial group representation theory; in particular, we summon the natural 
expansion in irreducible components and subspaces. Precisely, let V be any 
given vector space, on complex field C, and {Wg}gfzg a unitary representation 
of Q on V. Then V decomposes naturally as the Cartesian sum of irreducible 
subspaces V^: 

Be 



0vw 

,9=1 



0(dW®vW), (B.2) 



where the Vt'^J are the smallest invariant subspaces of V under the action of 



Ug. In (B.2), c is a scalar integer labeling which of the various irreducible 
representations (irreps) of Q is related to the given subspace y^^\ which, 
in turn, determines the dimension rhc of V^^^ Index c is commonly called 
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charge, or sector; for comfort, we will always label with c = the trivial 
representation vj"' = 1, of dimension rfiQ = 1, which always exists, regardless 
of Q. Moreover, in a generic representation, the same c-charged irrep can 
appear an arbitrary number of times Be] this (integer positive) Be number 
actually defines the degeneracy of sector c within representation Wg. 

In conclusion we can build a complete basis for V, in accordance to ex- 



pansion (B.2), listed as |c, mc), where G {l-.dc} labels degeneracy-space 



vectors, while rric G {l..mc} labels irrep-space vectors. Then the action of 
Wg is respectful of this decomposition: 



Wg\(p) = Wg i^Y^ (Pc,d.,mjc,dc,me) 

\ c dc=l mc=l / 
dc mc / fhc 



(B,3) 



c dc=l m'=l \mc=l 



or more simply, by exploiting the direct sum formalism of (B.2) we can write 

c 

\c] 

where 1 is the identity operator, and the unitary matrix Vg is the c-charge 
irrep of the group element g. An important remark is the following: undoubt- 



edly, the set of allowed c values in expansion (B.2) and (B.4) is a peculiar 



property of the representation W we chose, and similarly also the degener- 
acy numbers dc- Instead, the irrep dimensions rhc depend only the group Q 
itself (after irreps and charges c have been associated once and for all), not 
specifically on W. In particular, let us recall that if Q is an abelian group, 
then all its irreps have always dimension one, i.e. rhc = 1 regardless from c. 



B.2 Network geometry preservation 

We will now select an arbitrary network geometry, and focus on the mani- 
fold of Tensor Network states {|^tn)} which can be generated through that 
geometry. A preliminary argument that must be introduced, is that {I^E'tn)} 
is closed under the action of any pointwise symmetry group: 

|^') = f^ri^TN) =^ 1^') e {I^tn)}, (B.5) 

telling us that the target state |\E'') admits an exact tensor network repre- 
sentation (with the same graph design). This follows trivially from the fact 
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that the operation C/®^ is a tensor product of local and invertible transfor- 
mations. Therefore the target state |\&') must yield the same entanglement 
properties of the original state |\E'tn)- Another way to see the equivalence, 
from an algebraic viewpoint, is that every of the local Ug transformations 
can be adsorbed into the nearest linked tensor, thus effectively recovering 
the original network structure (and unaltered bondlink dimensions D). 



B.3 Zero charge Tensor Network states 

The first, and simplest, step to undertake if wc want match Tensor Network 
ansatzc and symmetries is to provide a characterization of TN-states which 
are invariant under the pointwise symmetry group, i.e. those I^^tn) which 
it holds 

= t^ri*™) e g. (B.6) 

It is clear that, in order for the unitary symmetry group to have no effect on 
them, these must strictly belong to the (9o-dimensioned) sector c — 

of Ug. Therefore, we will call states 1^^^) belonging to this sub-manifold, 
zero charge Tensor Network states. 

Clearly, if all tensors of |^xn) invariant under the action of Ug (applied 
at the physical links solely), the state I^'tn) invariant. This condition is 
sufficient, but not necessary. Indeed the application of Ug could map the 
tensors into a gauge-transformed version of the original ones, and the state 
would be unaltered anyway. In practice, tensors T of |^tn) J^ust be such 
that, for every g E Q, applying f/®^ is equivalent to performing a gauge 
transformation. For obvious reasons the selected gauges must form a group 
(intended as a subgroup of the whole gauge group), representing Q. But since 
gauge transformations of Tensor Networks are made of local, independent 
isomorphisms upon doubly-connected links, every link is implicitly carrying 
a group representation separately. 

In conclusion, we may associate to every non-physical link a a unitary 
representation W(^a) oi Q and a direction, and tensors T must be symmet- 
ric, i.e. invariant under the action of Q on all of their indices (with their 
respective representation and direction): 

/"phys \ r IV \ /"°"* \ 

rg*'= E nra: nK)j; n[w'(..)j:: ^t" 

{r},{v},{w} \ a / \ a ^" J \ a J 

(B.7) 

for every g E Q, where tensor T is connected to ctphys physical links, ctin 
incoming virtual links, and ctout outcoming virtual links. The representation 
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U is the only one fixed by the physics of the problem, while the W^(a) are 
arbitrary, chosen by the user. We are strongly encouraged to model the 
W^a) on the problem we are investigating, especially if we are implementing 



a simulation algorithm. It is instructive to rewrite equation (B.7) as follows: 





(B.8) 

The diagrammatic formulation of Tensor Network formalism will help us to 



understand why equation (B.7) is the standard requirement for a zero charge 



TN-state. Precisely, consider the following example: 






(B.9) 

where a three-sites state is involved, on which we are applying the pointwise 
transformation f/®^ for a given g E Q. The first equality in (B.9) is just 



creation of pairs operator-antioperator, which is freely allowed since it is a 



gauge transformation (as discussed in section [44). By choice, we are picking 



exactly the unitary operators W{^ot)g (green circle) and W^^^^ = W^^-^^ (orange 
circle), with the selected representation for each virtual link a. Now we 



apply the tensor symmetricity requirement of eq. (B.8 ), which automatically 



gives the second equality. Then the state is invariant under the action of the 
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whole pointwise symmetry group Ug^. We showed an example with a given 
network geometry, but it is obvious that this argument applies equivalently 
to every Tensor Network state whose tensors satisfy (B.7), (B.8). 

A peculiar ingredient in this framework is selecting the virtual represen- 
tations W(^a) of the symmetry group. It is true, as we stated, that for any 
choice of those representations the resulting TN-state would be zero charged; 
but at the same time, the description capabilities of the symmetric Tensor 
Network ansatz may depend on such choice, and typically will. 



B.4 Symmetric tensors fragmentation 

The prescription ( |B.7[ ) of using symmetric tensors in order to generate in- 
variant network states, can be interpreted [70] as a generalization of Schur's 
lemma to a linear operand (tensor) with an an arbitrary number of indices. 
Precisely, consider the case investigated by Schur, where an operator O com- 
mutes with a unitary representation 1^ of a compact symmetry group Q; the 
O can have nontrivial support only on the degeneracy space, i.e. 

[O,Wg] = =^ O = (OW ® l^^x^J , (B.IO) 



using the same irrep subspace decomposition of (B.2). It is easy to check 
that operators in this form are the only ones commuting with every Wg of 



equation (B.4). Then, as it was discussed in refs. [Ml [70] a similar argument 
applies to every symmetric tensor, whichever its amount of indices might be. 
Ultimately, symmetric tensors decompose in such a way that the degrees of 
freedom which are not fixed by symmetry (variational DOF) are isolated, 
and separated from the symmetry constraints (structural DOF). This de- 
composition takes place at every node of the network structure, substantially 
fragmenting a single tensor T into a fully-variational tensor R, and one (or 
more) structural-tensor S; thus actually splitting the original network into a 
pair of connected superimposed graphs. 

We will now explain how this symmetric tensor fragmentation scheme 
is performed in practice, by considering tensors with a limited correlation 
number (up to four attached links). We will proceed step-by-step starting 
from the simplest cases, and sketching diagrams whenever possible for clarity 
and comfort. 

Bondlink fragmentation - Before considering fragmentation of ten- 
sors, it is useful to understand how symmetry relations decompose network 
bondlinks themselves. Indeed, we mentioned that to every link a we as- 
sociated a representation W^a) of Q, as well as a direction which helps us 
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to discriminate between the application of the direct unitary transformation 
W(^a)g and its inverse W^^^^ (according to this picture, we may think open 
physical links as outgoing links). Then, every value j of the link a, intended 
as a virtual state | j)^ decomposes according to the irrep subspace expansion: 
\j)a \c,dc,mc)a- Therefore, the network link literally splits into: 



link a 



d, 



c 



^^^^^^^^^^m rric 

(B.ll) 

where the red fragment-link carries the charge index c, the green one holds 
the degeneracy index dc, and the blue one keeps track of the irrep basis 
vector label rric. As we stated, the number rhc of allowed values for the index 
rric depends on the value of c, in other words the blue link has a dimension 
depending on the value of the red link. Similarly, the degeneracy number 
9c"'* depends on both c and the representation chosen on link a. This is the 



purpose of the wavy arrows in (B.ll), to remember that the value of the red 
link influences the dimension of the blue and green links tied to it. The total 
bondlink dimension becomes then 

Dc.= Yl ^cy<di''\ (B.12) 

c e irreps(5) 

where we can sum over all irreducible representations of Q. with c being 
the charge. To control Da and keep it finite, even when there are infinite 
independent irreps for Q, it is sufficient to set the dc°'^ so that only a finite 
number of them are nonzero. 

To make a practical example, let us assume that Q is SO (3), then rhc = 
2c + 1 with c G N. Assume we are performing a renormalization process, 
keeping track of pointwise spin rotation symmetry; and three spin-1 sites 
were already renormalized into the given link a. Then, by spin-sum rules, 
we might want to describe for W(^a),g- one singlet, three triplets, two quin- 
tuplets, and one 7-plet, so that do = 1, di = 3, 83 = 2 and 84^ = 1 {Be = 
for c > 5). Clearly, if we take these degeneracies, we are keeping all the 
state information, and not actually renormalizing anything: indeed the total 
bondlink dimension D is 27, equal to d^. 
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As we stated, if the symmetry group Q we are considering is abelian, 
rhc is always 1, regardless from c or the link a; so the blue link- fragment 
allows only one value, and therefore is futile. Indeed in abelian symmetry 
frameworks, typically only the red and green sub-links appear, as they are 
the only ones needed. 

One-link tensors - a tensor having only one index behaves like a vector. 
This means that it is symmetric only if Q acts trivially on it. This implies 
that its support must be restricted to the c = sector, and j = Oq labels 
completely the space |0, 9o, 1), which is fully degenerate 

Two-links tensors - This is the case considered in Schur's lemma. 



Equation (B.IO) can be read in the following terms: a two-link symmetric 
tensor T, written in the basis \i){k\ — )■ |c, 5c, mc) (g, 9^, m'^l must preserve 
both charge and irrep label m. Therefore T naturally decomposes as 



Ti = X {SDZl with {S, 



c\mc 
q'm' 



^c,q ^m,m' 



(B.13) 



where R is the variational fragment of T, and S, the structural one, is equal 
to the identity 1. Let us sketch fragmentation diagrams for one-link and 
two-link tensors: 







(B.14) 

Notice that in the case of one-link symmetric tensor, we disregarded the irrep 
vector sub-link (blue one) since the only relevant charge is c = which is the 
trivial irrep, and thus 1-dimensioned. 



The peculiar double-delta form for S in (B.13) actually depends on the 



fact that we are considering one ingoing link and one outgoing, which is the 
most common setting. Other choices of directing links (say two incoming or 
two outcoming links) lead to different structure tensors: (5*^)^? = 6c,q ■ 
The best way to recover these setups is by contracting a 3-link symmetric 
tensor (we will discuss it shortly) with a 1-link one. 

The two-leg symmetric tensor case is the first where we encounter an 
actual reduction of variational parameters, by employing symmetries, while 
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Figure B.l: Fragmentation diagram for a three-link symmetric tensor T. 
The structure tensor S is made of Clebsh-Gordan coefficients for irreps of 
the group Q. We are assuming Q to be multiphcity free; otherwise, a third 
tensor fragment appears. 



preserving the original total bondlink dimensions D^, and thus entanglement 
features as well. Indeed consider respectively T and R of ( B.14[ b): 



T i — D^Dg descriptors, 



R 



^9].°lc}^^^ descriptors. 



(B.15) 

Then, since Da = Xlc^c^Wc, we typically end up with an effective amount of 
parameters -D^ ^ ^ D^Dp ~ way smaller than the original one, especially 

when the active sectors (c values for which Si"' > 0) are many. 

Three-links tensors - The tensor product of two irreps, with charges a 
and h respectively, is still a representation of Q] and thus can be decomposed 



in a direct sum of irreps according to (B.2) and (B.4): 



9c = l 



(B.16) 



where /x^^ is the number of copies of V^^^l appearing in the tensor product 
representation. 
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Here it is comfortable to assume that the group Q is multiphcity free, i.e. 
that fi'^ i^ can be either zero or 1, no matter the sectors. This is quite a typical 
case for physically relevant symmetries: SO (3) and SU(2) which take into 
account rotational invariance, as well as every abelian group that keeps track 
of particle number and parity conservation, are multiplicity free symmetries 
(but not SU(3), for example). In this framework, the Wigner-Eckart theorem 
provides a remarkable fragmentation law for a three-leg tensor T: 



where the structural tensor fragment S contains the Clebsh-Gordan coef- 



ficients for irreps of Q, which are well defined thanks to (B.16) and the 
multiplicity freedom requirement: 

('5a,(.):::L< = (c,mc|a,ml;6,m','). (B.18) 
An analogous decomposition with different structural tensors S holds for 



other direction configuration of the links connected to T. Figure B.l shows 



the diagram for (B.17), the case we considered. 

We want to check, in this three-link tensor scenario, the effective reduction 
of variational descriptors caused by enforcing symmetry relations. So, let us 
compare the amount of parameters: 

T^D^DpD,r^D\ R^D'''^J2Y1 E (B.19) 

a b c6{a©6} 

where the innermost sum spans the sole sectors c that are achievable by 
fusing together charges a and b, via (B.16). It is clear that D^^^ is very small 



compared to D^; it is more likely to scale like rather than D^, especially 
if the symmetry Q is abelian. 

Four-links tensors - The tensor product of three irreps VW, yt^] and 
yt'^l may contain several copies of the same irrep V'"^!, even when the symmetry 
group Q is multiplicity free. However, we can workaround this issue by fusing 
together two charges, say a and b, beforehand, operation which is well-defined 
through Clebsh-Gordan sum rules. 

Let us focus on a directing configuration having one ingoing link, and 
three outgoing ones. Let then e = a©6, where © stands for charge fusion, and 
clearly g = e © c. We have to take into account all the allowed intermediate 
charges, i.e. those e for which 



(B.20) 
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Figure B.2: Fragmentation diagram for a four-legs symmetric tensor T. 
The charges a and b fuse into a set of intermediate charges e, by means of 
Clebsh-Gordan tensor So- Index e also enters in the variational fragment R, 
to keep track of degeneracies arising from the product of three irreps. 



Adopting the intermediate charge scheme is sufficient to label and address 
separately the different copies of the same irrep V''']. Then the symmetric 
tensor T fragments as follows: 



rp i 



\9q 



q,e N 
b,cj 



m' ,m, ,m" 



(B.21) 



but, at the same time, we know that the structural fragment is obtained by 
two consecutive fusions, so in the end 



), (B.22) 



i.e. S further decomposes into a pair of sub-fragmented structure tensors 5*0, 
coinciding with the Clebsh-Gordan tensor we used for the three-links case, 
and we must contract over the m+ levels related to the e-charge irrep. The 



fragmentation scheme we just introduced is picted in figure B.2 



Of course, we could alternatively decompose T by fusing together b and c 
at a first step, and then fuse a with / = 6©c, as in figure B.3 b. In this case, 
we have different variational R and structural S tensor fragments, according 
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a) g , b) q 




c a 




Figure B.3: Two schemes to group together three irreps, of charge a, b, 
and c, by adopting the intermediate fusion charge framework, either e or / 
depending on which irrep is the last to fuse. The intermediate charge takes 
into account arising degeneracies in the irrep q. 



to: 



E 



E 



X r9^ 



(B.23) 



which is the right-left specular of the diagram |B.2 

The two intermediate fusion schemes we showed are in strict connection. 
Indeed, if we consider the partial fusion basis |e,m+) and \ f,mj) they are 
related [68] by the 6-index tensor F, (e.g. for the group Q = SU(2), F 
coincides with Wigner's 6-j symbols), so that: 



s: 



qj 



E« 



a,b,c\ 
d ) 



e y qq,e 



(B.24) 



At the same time, since the global tensor T is uniquely defined and does not 
depend on the choice of intermediate fusion basis, the variational fragments 
R and R must transform accordingly: 



*e 
f 



X R 



q,e 



(B.25) 



Theoretically, we could also study the scenario where a and c fuse together 
into an intermediate charge h, which fuses with b afterwards. But in order 
to investigate this framework, we should exchange the order of links, say a 
and b, before the first fusion takes place. As we discussed in section |4.3 



when swapping Tensor Network indices is needed, the inner nature of the 
degrees of freedom we are describing (i.e. whether they are spins, fermions. 
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bosons or even anyons) manifests as an exchange-statistic of links themselves. 
Therefore, the resulting fragmentation picture will depend on the statistic 
obeyed by the particles we are describing . 

Generalization of the equations we encountered in this section to other 
tensor topologies can be done by hand following the same intermediate charge 
rules. This concludes our discussion regarding Tensor Network fragmentation 
due to symmetry relations. 

B.5 Finite charge Tensor Network states 

The fragmentation schemes we just analyzed for different tensor topologies, 
derived from the argument that all the tensors in the network ought to be 
symmetric to make the global state I^'tn) invariant under the application of 
C/|^. Although this is an intriguing context, it is quite limited for practical 
purposes of studying physical settings. Indeed, in most variational problems, 
we are actually interested in working with an arbitrary, fixed, finite symmetry 
charge q. 

For instance, assume our model Hamiltonian H manifests a pointwise 
U(l) invariance, i.e. undergoes a particle conservation law: this means that 
global states, expanded in irreps subspaces as \c,dc,mc), belonging to dif- 
ferent sectors c are not coupled by TL. Then we might wish to achieve the 
ground state of restricted to a given charge, i.e. number of particles; this 
is definitely a physical question. Analogously, we could be interested in de- 
scribing the lowest energy state of a SU(2)-invariant T-L, on a spin-^ lattice, 
with total spin, say, |. 

Achieving a charge-selective ground state is no trivial task from a varia- 
tional point of view. It can not be performed blindly, by starting from the 
correct sector and hoping that the symmetry-invariant dynamics will prevent 
other sectors to be explored. This approach typically fails, as numerical er- 
rors will inevitably introduce other-sectors fiuctuations, which will ultimately 
break the symmetry, and push the algorithm towards the absolute ground 
state. The right way to deal with this charge-selective quantum problem is 
by forcing the variational wavefunctions to belong to the right sector. Here is 
were Tensor Networks succeed: we will be able to recover fragmentation rules 
we learned in the previous section, and apply them even for finite charge c 
TN-states, with c chosen by the user. 

Therefore, let |\E'tn) be our Tensor Network state, for a given graph ge- 
ometry, bondlink dimensions D, and consequent entanglement bounds. We 
will require that I^'tn), when expanded in irrep-subspace basis \c,dc,mc) of 
the pointwise symmetry group Ug^, has nonzero components only for a given 
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source 




^ TensorNetwork ^ 




SL 



Figure B.4: Diagram of the finite charge q Tensor Network prescription. 
The network is partially contracted into the violet dashed tensor, only the 
source node and the (structured) charge selector node are highlighted. 



charge q: 



dc rile 



Let us recall that dc are degeneracy indices, while are irrep basis indices. 
So when a symmetry transformation occurs, it interferes with the but 
leaves the dc unaltered, separately for every sector c. Then, our fixed-charge 
TN-state transforms as 



I^TN/ 



with vl''^ being the irrep unitary matrix, charge g, group element g. 

Now, we will give a prescription on the Tensor Network itself that will 
allow only states in the form I^I/tn), i.e. (B.26), to be variationally generated. 




[9] 



\q,dq,mq) 



(B.27) 



1. Choose a single tensor in the network, a node in the graph. We will 
refer to this tensor as source node. 

2. Direct the graph, i.e. associate a direction to every network link. Also, 
we require that every tensor in the network has at least one incoming 
link, except for the source node, which instead must have no incoming 
links. Such a directing scheme always exists, as long as the graph is 
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Figure B.5: The pointwise symmetry transformation f/®^ (pink boxes), 

when apphed to I^I/tn), hterally 'jumps over' the symmetric part of the net- 
work, which is formed by the original graph, including the source node. On 
the added link, it becomes Vg'^^ (green circle), as stated by equation (B.29). 



fully connected. Moreover, if there are no closed loops, this scheme is 
unique. Physical links are meant as outgoing. 



3. Associate a representation W(^a) of Q to every link 



a. 



Add a single one-leg tensor to the network {selector node), which con- 
nects to the source node. It is easy to see that TN-entanglement bounds 
are unaltered by this graph geometry adjustment. The link we just 
added is directed from the selector node to the source node. The rep- 
resentation of Q associated to the new link is the g-charge irrep Vg 



Symmetrize every tensor in the network, except for the selector node 
(but including the source tensor), according to directions and repre- 
sentations W(^a) of the links touched. This means that every tensor T 
decomposes into its variational R and structural 5* fragments, as we 
discussed in section [R4l 



6. Fix the charge selector tensor to be T^ 



c,mc 



Sc,qCmg, C arbitrary. 



This concludes the prescription. The Tensor Network state then reads as in 
figure B.4 The charge selector tensor can not be symmetric (unless Cm, = 0, 
or q = 0), as its only connected link has a nontrivial representation V^'^\ And 
in fact, I^I^tn) will not be Ug^ invariant since not every tensor is symmetric 
(actually, all but one are). 

Then, consider Z, the contraction of all the tensors in the final network 
except for the selector node. Clearly Z is symmetric, as contraction of sym- 
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metric tensors; i.e. it holds 

rfin d 



Kls, = E Y.^^i'^^U^rn, U^U,)s„r, ZT.l.r, ^9 G ^, (B.28) 

m'q {r) \ I / 

or, in matricial form, U®^ ■ Z ■ v!f^ = Z. In conclusion, when we apply the 
pointwise symmetry group to the Tensor Network state, we obtain 

d rhq 

t^r = uf' E E i^TZs.c^.) •• • sl) = 

S1...SL mq 

d rhq / \ 



S1...SL mq 



9 



which is formally equivalent to (B.27), thus proving that our prescription is 
sound. 

The reason why we needed to direct the graph so that it had a single 
source of directions (the source node) is to spread the information about q 
to the whole network. If this is not the case, then one can identify regions of 
the network insensitive to q thus actually behaving like symmetry-invariant 
zones: e.g. a party of sites which are always empty. Although the resulting 
state would still be a g-charge Tensor Network, it would be far more trivial. 

With the construction we just introduced, we are finally able to under- 
stand and exploit the strict relationship that ties symmetries and Tensor 
Networks through representation theory. In these sections we developed se- 
lection rules and manipulation techniques to embed symmetries into Tensor 
Network variational ansatze, allowing us to address charge-specific problems, 
and to meet a drastic speed-up in computational time. 



B.6 Example: Symmetries in MPS 

We would like to conclude this appendix chapter by applying the symme- 
try arguments and techniques upon a most common template in the family 
of Tensor Network, namely on Matrix Product States, the variational coun- 
terpart of DMRG algorithms. Methods for dealing with symmetries within 
the Density Matrix Renormalization Group framework were already known 
before the acknowledgement of Tensor Network states |5] , still, the in-depth 
understanding of both MPS representations and symmetric Tensor Network 
states, allows us to build a formulation for finite charge MPS which is com- 
pact, elegant, and efficient. 
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Here we will work with open boundary conditions MPS, as the no-closed 
loop geometry encounters less accidents, and show the MPS-fragmentation 
scheme respectively for an abelian symmetry group Q, and then for a non- 
abelian one. Generalization to PBC is not trivial but possible nevertheless. 



B.6.1 MPS with pointwise U(l) 



The abelian symmetry U(l) has infinite one-dimensional (m^ = 1, Vc) non- 
equivalent representations labeled by integer numbers c G Z. U(l) is used to 
take care of particle conservation, when the Hamiltonian "H has only terms 
that preserve particle number; indeed its fusion rule © corresponds to the 
simple sum of two integer numbers, i.e. c® c' = c + c'. 

A natural way to choose the source node (defined in the previous section), 
in order to characterize MPS states with an arbitrary particle number q, 
is to choose one of the edge blocks, say the one at right boundary. MPS 
tensors have three connection links, so we can use three-leg fragmentation 



rule (B.17) to split a block into structural and variational part. Also, recall 
that since rhc = 1, we have no need for blue (irrep vector) links. Ultimately, 
the resulting fragmented-MPS reads: 





(B.30) 

where the structural tensors are homogeneously defined 5'cj_i,cj = ,c _i+s • 
The yellow tensor Q is the charge selector node, properly connected to the 
source node A^^'^; its purpose is to select the global sector: = ^c^,? with a 
total charge q chosen by the user. The tensor fragments i?'-'' are completely 
variational, and we can freely manipulate their parameters, for instance, to 
lower the total energy, without constraints: we will always remain forcefully 
in the correct q sector due to the presence of 5* fragments. 

An intriguing feature of this abelian symmetric-MPS is that it is always 
operationally possible to gauge-transform it into the left (or right) gauge. 



while preserving the fragmentation scheme (B.30). The basic idea is to per 
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form a singular value decomposition of R'^^ tensors separately for every c,-: 



i^}) r.y-n = E = A,., (Kj'lt)^^ (B.31) 

i-1 



-1,9. 



/3c 



This is equivalent to performing an SVD of a block diagonal matrix, by 
actually singular value decomposing every diagonal block separately. This 
is not only formally meaningful, but also cheaper in terms of computational 
time. In the end we can recover all the engineering we developed in section 



2.9 , which strongly exploited left and right gauges, and further enhance its 
computational power by embedding symmetries. 



B.6.2 MPS with pointwise SU(2) 

The Heisenberg model a, ■ aj is the archetype of an SU(2)-invariant lattice 
Hamiltonian. Being a continuous symmetry, SU(2) can not be spontaneously 
broken, so its ground state has to be a total spin (provided it is possible 
by fusion rules). Nevertheless, we might be interested to describe either this 
ground state, or maybe the lowest energy level at fixed total spin q. 

SU(2) irrep charges c are typically labeled by integer and half-integer 
positive numbers, i.e. c G lNl/2, and corresponding irrep (blue link) dimension 
rhc = 2c + 1. The MPS fragmentation scheme [71] then reads: 
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(B.32) 

where structure fragments S are Clebsh-Gordan coefficients. The tensor Y 
can be any random tensor; no algorithm based on a SU(2) invariant bench- 
mark can determine or variate Y because it is the only non SU(2) invariant 
component of the MPS network. 
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